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Prologue 



Motivation 

In astrophysics a multitude of systems and configurations are described with concepts 
from hydrodynamics coupled with gravitation, radiation and magnetism. Mathemati- 
cally radiation hydrodynamics (RHD) and magnetohydrodynamics (MHD) are systems 
of coupled nonlinear partial differential equations. Multiple fields in astrophysics have 
been adopting and developing sophisticated numerical methods for solving these PDEs 
in various applications. 

Explicit numerical schemes in computational fluid dynamics have been popular for 
many years and received substantial boosts due to advances in technology and paral- 
lelizing techniques. Riemann solvers and related methods are extensively adapted in 2D 
and 3D computations however they have one main disadvantage. An inherent limitation 
for time steps in explicit schemes impedes applications where various time scales are of 
interest. Hence a majority of these calculations emphasize on small temporal and spatial 
scales but fail to treat phenomenons properly on very diverse scales. 

Techniques for implicit numerics are computationally more expensive since they are 
hardly parallelizable but do not limit the increments of the scheme intrinsically. Implicit 
RHD in ID has proven favorable for describing astrophysical processes on large spatial 
and temporal scales like nonlinear pulsation or long term development of supernovas. 
The idea of this thesis was to study suitable generalizations of these concepts to 2D 
and 3D. Additional degrees of freedom in multi-dimensional implicit RHD would dis- 
close the treatment of large scale convection and how it interacts with rotation, mixing, 
the coupling of rotation and pulsation, transport of angular momentum and mass loss, 
interactions of discs and stars during star formation processes etc. The target of this 
thesis was to find the strong conservation form of the equations of radiation hydrody- 
namics in non-steady curvilinear coordinates and to study multi-dimensional adaptive 
grid generation. 
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Synopsis 



The Euler equations of hydrodynamics, the Maxwell equations as well as radiative trans- 
port equations are hyperbolic PDEs that connect certain densities and fluxes via con- 
servation laws. Generally they emerge from natural symmetries that constitute major 
principles in mathematical physics. The numerical implementation of these equations 
essentially needs to comprise these qualities. Applied mathematics have developed an 
articulate framework of numerical methods for conservation laws that ensure the con- 
servation of mass, momentum, energy etc. if applied properly. The main challenge is 
to compute the fluxes correctly which indeed will be the issue of the first three chap- 
ters in this paper. Self gravitation is described via the Poisson equation that poses an 
extra challenge due to its elliptic nature. While ID computations avoid solving this 
pure boundary value problem by considering an integrated form we will have to design 
a different approach to self gravitation in 2D and 3D. 

A generalization of implicit conservative numerics to multiple dimensions requires ad- 
vanced concepts of tensor analysis and differential geometry and hence a more thorough 
dedication to mathematical fundamentals than maybe expected at first glance. Hence 
we begin to discuss fundamental mathematics and physics of RHD with special focus on 
differential geometric consistency and study numerical methods for nonlinear conserva- 
tion laws to gain a solid definition of the term conservative. The efforts in tensor analysis 
will be needed when applying Vinokurs theorem to gain the strong conservation form 
for conservation laws in general curvilinear coordinates. Moreover, it will be required to 
slightly reformulate the artificial viscosity for such nonlinear coordinates. 

Astronomical objects are characterized by fast flows and high propagation speeds on 
the one hand but astronomical length and time scales on the other hand. Implicit nu- 
merical schemes are not affected by the Courant Friedrichs Levy condition which limits 
explicit schemes to rather impracticably small time steps. Implicit methods however pro- 
duce algebraic problems that require matrix inversion which is computationally expen- 
sive. In order to achieve viable resolution, adaptive grid techniques have been developed. 
It is desired to treat processes on small length scales like shocks and ionization fronts as 
well as physics at the extent of the objects dimension itself like large scale convection 
flows and pulsations. The combination of implicit schemes and adaptive grids allows to 
resolve astrophysics appropriately at various scales. 

In the last chapter of this paper we study problem oriented adaptive grid generation 
in 2D and 3D. We establish three main postulations for an ideal grid and analyze several 
feasible approaches. 
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1 Fundamental Mathematics of RHD 



The equations of radiation hydrodynamics (RHD) are a system of nonhnear second 
order partial differential equations with hyperbolic and elliptic parts. Generally they 
are posed as an initial value problem with boundary conditions. The simultaneous 
solution of hyperbolic terms, as coming from the dynamics and the strictly elliptical 
Poisson equation constitute one major problem for finding suitable solution methods. 
In the following chapter we will retrieve the fundamental mathematics of hyperbolic 
conservation laws as they appear in hydrodynamics and radiation transfer as well as 
some basic concepts for elliptic partial differential equations like the Poisson equation in 
order to motivate certain numerical concepts that will arise in chapter 3. As we will see, 
accurate analysis provides expedient approaches to our numerical problem. 

1.1 The Hyperbolic Part 

Dynamical physical processes are mathematically described via hyperbolic partial dif- 
ferential equations^. The wave equation, the Burgers' equation, the equations of hy- 
drodynamics and the Dirac equation are prominent examples of hyperbolic evolution 
equations. They always characterize systems with finite propagation speeds, a quality 
often referred to as causality. Information from a point in space and time is causally 
associated solely with a finite region, a domain of dependence. In chapter 3 we will resort 
to this important issue when we discuss stability and convergence of numerical methods. 

In the following section we want to define some fundamental mathematical terms and 
present important concepts that will accompany the entire paper. 

1.1.1 Hyperbolic Conservation Laws 

The Euler equations of hydrodynamics as well as the angular moment equations of 
radiative transfer can be derived from conservation laws originating from one of the 
most fundamental theorems of mathematical physics, the Noether Theorem (see e.g. 
Bhutani (1981) [ I], Kambe (2007) [21]). In hydrodynamics, the conservation of mass, 
momentum and energy leads to the continuity equation, the equation of motion and 
the energy equation. Analog correspondences are derived for radiation. Mathematically 

^The Sdirodinger equation of quantum mechanics is excepted. 
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1 Fundamental Mathematics of RHD 



the equations of radiation hydrodynamics form a system of hyperbolic conservation laws 
that describe the interaction of a density function d 



d : 



and its flux f 



f : U 



X [0, oo) 



U C 



U open. 



The temporal and spatial change of the integrated density in a connected space S 
then equals the flux over the boundary dQ, where S is the outward oriented surface 
element. 

dtJddV + Jf-dS = Vt>0 (1.1) 
n an 

Hyperbolicity is granted if the Jacobian matrix associated with the fluxes Vjf has real 
eigenvalues and there exists a complete set of eigenvectors. Provided it is physically jus- 
tified, demanding this quality on a set of equations even can help to formulate constraints 
to the solutions (see subsection 2.2.3). 

For the in section 2.1 further motivated physical problem of hydrodynamics this system 



d(x,t) 



an example^. 






/ P(x,t) \ 




/ pu 


(pu)(x,t) 


f(d) = 


upu + P 


V (pe)(x,t) / 




\ peu + P • u 



(1.2) 



dtj pdV + J pu-dS 

n an 

dtJupdV + J {upu + P)-dS 

n an 

dtJpedV + J {peu + P • u) • dS 



(1.3) 



Let f be a continuously differentiable function, then we can rewrite equation (1.1) via 
the divergence theorem as follows. 

{dtd + dwf{d))dVdt = Vt>0,17cM" (1.4) 

t n 



^Notation: Summations respectively projections of tensorial quantities are always denoted by • or : 
whereas tensor products (g) are omitted, i.e. uu = u ® u. 
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1.1 The Hyperbolic Part 



Locally the system of partial differential equations 

atd + divf(d) = (1.5) 

gives pointwise solutions for our vector field d that describes the conserved densities or 
state variables of mass, momentum and energy. With initial condition (or the initial 
state of the density function) do 

d(x,0)=do(x) VxER'^ (1.6) 

Equation (1.5) with initial condition (1.6) is called Cauchy problem. The main theorem 
for this class of partial differential equations associated with analytic initial conditions 
is the theorem of Cauchy-Kovalevskaja which states local existence and uniqueness also 
for nonlinear Cauchy problems. However, with analyticity demanded this theorem has 
hardly any practical relevance. Questions of local and global well-posedness for general 
nonlinear Cauchy problems require sophisticated concepts of functional analysis (see e.g. 
Adams [i], Evans [14]). We want to mention some basic concepts of these variational 
formulations of partial differential equations in an extent that they arise in the context 
of radiation hydrodynamics and accordingly numerics. 



1.1.2 Weak Solutions 



Since even the simplest examples of one-dimensional scalar conservation laws like the 
Burgers' equation have classical solutions only in some special cases, one has to broaden 
the considered function space of possible solutions (the adequate function spaces are 
Sobolev spaces, see Evans [14]). For so called weak solutions, we appeal to generalized 
functions where also discontinuities can be dealt with. The generalized concept of differ- 
entiation of distributions shifts operations to test functions 7(x, t). These test functions 
7 : G G M" X M+ — M have compact support, meaning that there exists a compact 
subset (here i.e. closed and bounded) K such that 7(x, t) = 0\/x £ G \ K{j). 

(atd + divf(u))7dydt = V7GG (1.7) 

t>OR" 
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1 Fundamental Mathematics of RHD 



To simplify matters, we let these test functions be continuously differentiable. Then 
with partial integration and our postulations for 7 



dtd jdt = 7 d 



d dtj dt 



00 

7(x,o) do(x) - J ddadt 



divf(d)7(iy = 7f 



f grad 7 dV 



(li 



we derive the weak formulation of our conservation law (1.4). 



y" y" (d da + f (d) grad 7) dVdt = - J 7(x, 0) do(x) dV (1.9) 



t>OR" 



The function d G C°° is called weak solution of the Cauchy problem (1.4), if it satisfies 
(1.9) and d £ U with dg G C°° . This weak solution is not necessarily unique and 
normally further constraints have to be imposed. 



1.1.3 Rankine-Hugoniot and Entropy Conditions 

For the most physical problems it will be sufficient to look for weak solutions in the 
function space of piecewise continuously differentiable functions. The physical variables 
d are weak solutions, if they are a classical solutions wherever they are and at discon- 
tinuities (shocks) they need to satisfy additional conditions. As we will see in chapter 2, 
these conditions posses direct physical relevance when applied to a physical system like 
radiation hydrodynamics. The full derivation of the one-dimensional Riemann-prohlem 
for (1.4) considered can be found in the Appendix 6.1.2. 

We conclude, 

00 00 00 

[dda + f{d) dx^) dxdt = — J j{x, 0) do{x) dx + 

—00 —00 
00 

+ J j{x,x/ns) (^ ^(^^)-^(^-) _ (d, -4)) dx (1.10) 
* 

is weak solution of the Riemann-problem, if 

m) - fjdr) _ /^^ _ ^ X ^ Q 
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1.1 The Hyperbolic Part 



respectively 

fjdl) - fjdr) 

We call Us shock velocity and (1-12) Rankine-Hugoniot condition. 

As already mentioned, weak solutions of the Cauchy problem are not necessarily unique 
and thus may be physically pointless. However, there is a way to pick physically valuable 
solutions out of multiple. One method is to implement an artificial viscosity term 

^td* +divf(d*) =ei/Ad*, e>0 (1.13) 

and evaluate the equation for e — t- as limit case of the original one. The idea is 
motivated from physical diffusion, which broadens sincere discontinuities to differentiable 
steep gradients. The physical solution of the weakly formulated problem has to be the 
zero diffusion limit of the diffusive problem. However, in analytical practice this limit 
is rather bulky to calculate, hence simpler conditions have to be found. The common 
technique to do this is motivated from continuum physics as well, where an additional 
conservation law holds for the entropy of the fluid flow as long as the solutions remain 
smooth. Moreover, it is known, that along admissible shocks this physical variable never 
decreases, so the conservation law for the entropy can be reformulated as inequality. 

We regard a scalar entropy function o'(d) and an entropy flux i;^(d) which satisfy 

(9to-(d) + div (/>(d) = 0. (1.14) 

Assuming differentiable functions, we rewrite that conservation law via the chain rule 
and compare it to (1.5) multiplied with to obtain 

Vdaatd + Vd^Vxd = 
VdcrSfd + VdO-Vdf Vxd = 

Vdcxdivf = div(l) (1.15) 

where in the more dimensional case the appearing matrices of gradients have to fulfill 
some further constraints (see e.g. Godlewski and Raviart (1992) ['■']). For scalar equa- 
tions, it is always possible to find an entropy function of that kind. Furthermore we 
postulate convexity for the entropy function.'^ 

V^o->0, Vd (1.16) 



^In physical entropy one would require an eventual inequality with > in place of < but mathematical 
literature commonly chooses a" > 0. 
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1 Fundamental Mathematics of RHD 



To get our actual entropy condition, we rewrite (1-14) in the viscous form 

dta{d*) + div (/.(d*) = euVd<T{d*)Ad* (1.17) 
and integrate over an arbitrary time interval [to,ii]- 

J J (5t(T(d*) + div(/.(d*)) dVdt = J J VdCT(d*)Ad* dVdt = 

to Q to Q 

tl tl 

= j j (Vxd* Vdcj(d*)) • dSdt-eu J J V^.,d* VV(d*) V^^id* dVdt 

to dfl to C >o 

(1.18) 

When we now consider our non diffusive limit for e — 0, the first term on the right hand 
side vanishes without further restriction whereas the second term has to remain non 
positive. With partial integration and divergence theorem we get our entropy condition 
(1.19). 

tl 

J J (^dta {d) + div (j){d)'^dVdt < 

tl 

J a{d{x,ti)) dV < J a{d{x,to)) dV - J J ^{d) dSdt 



to n 



to an 



[1.19) 



For bounded, continuous pointwise solutions d* with d* — > d for e — t- 0, the vanishing 
viscosity solution d is weak solution of the initial value problem (1.4) and fulfills entropy 
condition (1.19). Generally spoken, applying the entropy condition systems with shock 
solutions unveils those propagation velocities that ensure that no characteristics rise from 
discontinuities which would be non-physical. For details and proofs we refer to LeVeque 
[26]. 

In chapter 3 we will discuss these analytic concepts in numerical context and work out 
its concrete applications for (radiation) hydrodynamics. 



1.2 The Elliptic Part 

As already elaborated radiation hydrodynamics are partial differential equations of the 
hyperbolic type. However, the dynamics are only one part of the challenge emerging in 
astrophysical applications. Self gravitation in classical astronomical objects is described 
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1.2 The Elliptic Part 



via the Poisson equation, an elliptical partial differential equation which emerges as 
pure boundary value problem for the gravitational potential. Unlike hyperbolic PDEs 
elliptical problems mostly describe equilibrium physics and their solutions are typically 
quite smooth as we will see. 



1.2.1 The Poisson Equation 

The Poisson equation describes the spatial distribution of a scalar field ^ 

$ : [/ C M" ^ M (1.20) 

induced by a source p 

p:UOW^R (1.21) 

and has the following form; the choice of the minus sign originates from the customary 
definition for ellipticity of differential operators. 

-A$ = p. (1.22) 

When it comes to our physical application of self-gravitation, we will pull the sign into 
the definition of the gravitational potential and consider A^g^av = ^T^Gpm- 

The Poisson equation is posed with Dirichlet or Neumann boundary conditions which 
means either boundary values or their derivatives are given. Often it is useful to impose 
mixed boundary conditions, where the boundary of the region is separated in a part 
with Dirichlet conditions 9$!^"^ and a part with Neumann boundary conditions d^^™ 

and dn = dn^'' u dn^"''. 

n-V^ = V^"" in5Q^^" (1.23) 

The Poisson equation can be understood as the inhomogeneous generalization of the 
Laplace equation A<I> = 0. Physically, latter equation describes a solenoidal, irrotational 
vector field u = V<I> with divu = rotu = 0. This circumstance implies that $ must 
be invariant under rotation and can be used to derive the fundamental solution. With 
Green's function and convolution of the inhomogeneous term p we derive the analytical 
solution of the Laplace equation^. In the light of these presumptions, we can investigate 
fundamental attributes of solutions of the Laplace equation in order to learn something 
about the solutions of the Poisson equation. 

*"I>(x) = $(r), diT = ^ . . . — > Green's functions of the fundamental solution of the Laplace equation 
in 3D is A(3) {-^) = S^'^^x) and $ = / G(x,x')p(x') dV is solution of A$ = p. 
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1 Fundamental Mathematics of RHD 



1.2.2 Harmonic Functions and the Maximum Principle 

One important aspect of the Laplace operator as it emerges in Poisson and Laplace 
equation is the fact that it acts regularizing; the according theorem follows. Hence it 
will be adequate to deal with differentiable functions when we look for solutions of (L22). 
Any function satisfying A<t = in [/ and is called harmonic function. 

Given $ is harmonic, we can derive a mean-value formula which states that on a 
spherical subset'^ i3(x, r) C U, a solution ^>(x) equals both the average of $ over the 
sphere dB and the average of ^ on the entire ball B. 

^(^^) = J^±}1^ f |,d5 = Eit+ll f ^dV (L24) 

dB{x,r) B{x,r) 

In 3 dimensions this simply yields ^(x) = J ^ dS = J ^ dV with correspond- 
ing integration ranges. In the light of these conclusions, we deduce some interesting 
properties of harmonic functions. 

The (strong) maximum principle states for a $ G C'^{U) nC{U) harmonic within U, 
i.) max^ ^ = maxg(7 <I>. 

a.) Furthermore, if U is connected and 3 xq S [/ such that u{xq) = max^/ it follows 
that $ is constant within all U . 

Assertion i.) implies that the maximum will be found at the boundary of the domain and 
a.) states that if a maximum is on the inside of the boundary the function $ consequently 
is constant within U . Another important conclusion of this theorem concerns the Poisson 
equation. It can be shown that for the boundary value problem 

— A$ = p InU 

^ = g indU (1.25) 

there exists at most one solution ^> G C2(f/) n C{U). With this uniqueness theorem we 
gain at least some comfort for the gravitational part of our numerical problem. As a 
further side note we mention the regularity theorem. It reveals that if <I> E is harmonic 
then it automatically is infinitely differentiable, i.e. ^> S C°°. 

Without dwelling on the Laplace equation any further, we recapitulate that math- 
ematical analysis of the elliptical part of our system of equations brings forward less 
sophisticated problems than the hyperbolic part. Obviously, uniqueness and regularity 
are soothing qualities when it comes to solving partial differential equations. All the- 
orems and proofs of this section can be found in Lawrence C. Evans (1998), p. 17 ff. 
[14]. 



^Closed n-dimensional ball -B(x, r) with center x, radius r > 0. 
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2 Fundamental Physics of RHD 



Astronomical objects are basically self-gravitating aggregations of (gaseous) material 
with magnetic fields and radiation. Radiation hydrodynamics as we will discuss them, 
factor out magnetic effects and solely consider radiative contributions through photons. 

The fundamental equations of hydrodynamics describe the motion of fluids in in terms 
of continuous quantities. The dynamics and kinetics of fluid particles are described by 
macroscopic functions like the mass density p, its velocity u, internal energy and gaseous 
pressure as primary variables mostly. In this model, infinitesimal volume elements of the 
fluid are still big enough that intermolecular forces can be neglected, meaning that they 
contain a sufficient number of molecules. The thermodynamic variables intrinsically 
are assumed to be well defined (equilibrium thermodynamics) and characterized by an 
equation of state. 

The radiative transfer equation describes energy transfer in form of electromagnetic 
radiation, affected by absorption, emission and scattering. The fundamental quantity 
to describe the radiative field is the spectral intensity /-y, which can be understood as 
macroscopic quantity to describe the photonic energy fiow. The mass density p is also 
source for the gravitational potential In the following sections we will physically 
motivate the set of equations of radiation hydrodynamics with gravitation bit by bit. 

One major intention in the upcoming paragraphs will be rigorous and consistent (in 
terms of differential geometry) definitions of the physical variables which will be crucial 
for the numerical methods in general cuvilinear coordinates investigated in section 3.2. 



2.1 Hydrodynamics 

As anticipated, the equations of hydrodynamics can be derived from conservation laws 
of mass, momentum and energy. Since the fundamentals of hydrodynamics are more 
than well documented (e.g. in Landau, Lifschitz [:. .]) we confine our portrayal to a brief 
synopsis and concentrate on a few interesting details. 
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2 Fundamental Physics of RHD 



2.1.1 Continuity Equation 

The conservation of mass is physically justified as long as no mass is generated or anni- 
hilated 

J pdV = (2.1) 



dm d 
dt dt 



V(t) 

and leads to the continuity equation. 



I {^dtp + div (pu)') dV = 

I dtp+ J {pu)-dS = (2.2) 



V{t) dV{t) 

Clearly the domain Q of the conservation law, introduced in section 1 . 1 is an n-dimensional 
Volume that is a bounded subset of the M", n = 1, 2, 3. The mass density p is a non neg- 
ative scalar field p(x, t) : M""*"^ — )■ M"*" and the velocity u a vector field u(x, t) : M""''^ — )■ 
M"'*'^. Solutions of these primary variables p, u, that is densities of the conservation 
laws, will emerge as piecewise continuous functions. 

As a side note the flow is called stationary if dtp = and incompressible if divu = 0. 
However, in astrophysical applications these special cases will fail to describe the physics 
appropriately. 



2.1.2 Equation of Motion 

In order to yield the equation of motion from the conservation of momentum it is re- 
quired to specify which forces act upon the fluid. Now we are concerned with a non-scalar 
conservation law and as we will see in chapter 3, in this case we will have to avail our- 
selves of general tensorial notation. Thus, at this point we introduce one main tensorial 
quantity (and secondary physical variable), the non-relativistic viscous stress tensor. 

T = puu -h P - cj 

= pu'u^ + P'^ - a'^ (2.3) 

T(x,t) : M"+i X ]R"+i M"+i x M"+i. This symmetric^ tensor {T^i = T^') describes 
contributions from momentum pu^ in j-direction, generally anisotropic pressure P^^ and 
friction o"*-^ forces. Diagonal elements express pressure and off-diagonal elements shear 

is symmetric, if conservation of angular momentum px x u is invoked, see derivation in the Appendix 
6.1.1. 
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2.1 Hydrodynamics 

stresses. The equation of motion in then yields^ 

J (9t(pu) +divT) = 

vit) 

J dt{pu)dV+ J (puu + P - a) • dS = 0. (2.4) 

V(t) dV{t) 

Special cases for the energy-stress tensor in fluid dynamics are e.g. the dust model where 
no particle interaction is assumed P = o" = 0. Anyway, in numerous physical applications 
the pressure forces can be expected to be isotropic. If so, the pressure tensor simplifies 
to a diagonal tensor which consists of the scalar gaseous pressure p and a measure of 
uniformity, the the metric tensor g, which yields P = pg. The divergence of that tensor 
can be further simplified to the well known gradient of the gaseous pressure. We will 
extensively occupy ourselves with the viscous pressure tensor a in succession, especially 
in section 3.2.4. At this point it shall only be mentioned that we neglect physical viscosity 
since the gas of a star is easily approximated as ideal fluid. As indicated in section 1.1.3, 
we will impose artificial dynamic viscosity when we concern ourselves with mathematical 
methods for solving nonlinear conservation laws. 

Implicitly we already have defined a thermodynamic equation of state which describes 
the state of the fluid in terms of pressure p, volume V and temperature T by imposing 
an ideal fluid. Since stars are largely composed of hot, low-pressure hydrogen gas, it is 
quite reasonable to apply this model. 



2.1.3 Energy Equation 

The energy balance of a fluid includes kinetic and pressure parts as well as inner energy. 
Latter is once more a thermodynamic quantity which is associated with the equation 
of state. We deflne a new non-negative scalar field, the specific inner energy e(x, t) : 
]gn+i _^ -^j^jc]^ jj^ CQQQ of an ideal fluid equals its thermic energy. The conservation 
of the total energy of a bounded fluid element in this strictly hydrodynamic case yields'^ 

j {dt{]^pv? + pe) ■M + peu))dV = 

V{t) 

j dt{]^pv? + pe)dV+ j ((^/9u2 + ^e)u+(P-a) -u) -dS = 0. 

v{t) av{t) 

(2.5) 

^In general relativity the conservation of momentum can be found in its most elegant form V • T = 0. 

In this case all variables are defined in a four-dimensional spacetime, where the 0-component of the 

covariant nabla operator contains the temporal derivative. 
^Using V[\v?) = uV • u. 
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This total energy equation can be dissected into a mechanical and an internal part, where 
former describes the change of kinetic energy of a fluid as the rate of work done on an 
element by pressure and stress forces. The gas energy part can be formulated as entropy 
conservation {dt{ps) + div (psu) = 0) and hence implies, that the flow is adiabatic. 

2.2 Radiation 

The fundamentals of radiative transfer and its moments can be found in several appro- 
priate text books (see e.g. Mihalas and Mihalas [32]), therefore we will concentrate on 
some noteworthy highlights in this section as well. In the previous paragraph we treated 
the physics of fluids which we now want to extend to radiating fluids in a nutshell. 

Radiation contributes to the total energy, momentum, stress and pressure in the fluid. 
In succession we will work out means of coupling among these physical quantities in the 
material. When we concerned ourselves with the fundamental mathematics of conser- 
vation laws, we denoted that also the fundamental dynamics of radiation will appear 
as hyperbolic conservation laws. In fact, the crucial equation to describe radiation is 
the radiation transfer equation (RTE), which basically expresses conservation of energy. 
However, we shall begin to define the fundamental dynamical properties of the radiation 
field. 

2.2.1 The Radiation Field and its Momenta 

The macroscopic quantities of the radiation field are based on statistical physics of a 
microscopic description of ultra-relativistic massless particles, i.e. photons. To establish 
the nexus to the continuum view which we will pursuit, we define the scalar non- negative 
photon distribution function /^(x, t; which can be interpreted 

as photon density in the corresponding phase space^. The central physical quantity for 
the upcoming description of the radiation field is the spectral intensity /-^(x, t; n, z^) 
defined by 

A = -^fT (2-6) 
That implicitly describes a radiative energy transport in the following sense 

d£^ = I-ydS ■ udojdvdt. (2-7) 

The spectral intensity at position (x,t), traveling in direction n with frequency v is 
defined to account for the energy 8^ transported by radiation of frequencies dv across a 

*The momentum of a photon, traveling in direction n is given by p-, = ^n. v is the frequency, c the 
speed of light and h the Planck constant. is invariant under Lorentz transformations. The subset 
7 is chosen to discern mere radiation (here photonic) variables. 



20 



2.2 Radiation 



surface element dS in time dt into a solid angle du around n. Further useful properties 
are the angular moments of the radiation field, that is the mean intensity (zeroth 
moment), the radiation flux F-,, (first moment) and the evidently symmetric radiation 
pressure tensor (second moment). The velocity of our particles will simply be given 
by the speed of light and their direction of motion u = cn, which also illustrates the link 
to the procedure we carried out for deriving the equations of hydrodynamics. 



J^(x, t;i/) = — I I^du) 



F^(x, t]v) = j I^n 



dijj 



P^(x,t;i/) = - I I^nnduj (2.8) 



c 



Analogously to our fundamental equations of hydrodynamics, we can derive dynamic 
equations for the radiation field by considering the conservation of density and mo- 



mentum Fy. 



J (dtJy + dw{cFy)^dV = 

vit) 

J (^dtFy + div {c^Py)^ dV = (2.9) 

V{t) 

From the viewpoint we pursued in section 1 . 1 discussing hyperbolic conservation laws we 
define a density function d(x, t;i/) = (J^,F^) and a fiux f(d) = c(F^,cP^). Equations 
(2.9) can hence be interpreted as homogeneous outlook to the moments of the radiation 
transfer equation with collision terms (2.12) in the forthcoming section. 



2.2.2 Interaction with Matter 

For the treatment of radiation hydrodynamics it will be necessary to describe the inter- 
action of the radiation field with the fiuid. Again the macroscopic approach chosen here 
to describe this interaction can be understood canonically with concepts from statistical 
physics and quantum mechanics. Moreover, it will be required to process the variables 
of the radiation field into the fiuid frame of the material correctly via Lorentz transfor- 
mation and account for the frequency change with respect to the lab frame that comes 
with Doppler shift. 

When radiation passes through material, energy generally will be removed depended 
on a material property, the opacity or extinction coefficient Xy On the other hand, the 
material releases radiation in terms of an emission coefficient t]^ which will be isotropic 
in its rest frame. Both properties are defined as non- negative scalar functions of posi- 
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tion, time, radiation direction and frequency (^7, f?7)(x, t; n, z^) : M^+^ x M'^ x M — > M+. 
Both processes, extinction and emission, can be further divided in thermal and scatter- 
ing events. Thermal extinction and emission converts radiation energy to thermal and 
vice versa, whereas scattering in principle changes direction and frequency of traveling 
photons. The amount of energy removed along a material element of length dl is 

d£^ = (x^^""^"^ + xT^)!^ dSdldudvdt (2.10) 

and the amount of radiation energy released by emission is 

d£^ = {rf^"'"^ + r?^^^*) dSdldujdudt. (2.11) 

The conservation of spectral intensity (2.9) considering emission and extinction yields 

J (^dfl^ + div (Iju) - r]^ + x-ylj) dV = (2.12) 
V{t) 

and is called radiation transfer equation in the laboratory frame. In practice of RHD, 
upper equation and its momenta have yet to be transformed into the fluid frame as 
mentioned and simplified in several ways. 

Firstly the monochromatic quantities I^,x-y, - ■ ■ will solely be considered with their 
collective effects by integrating them over all frequencies^. 

c - 1 f 
F^(x, t) = 47rH.y(x, t)= J L^ndudv 

P^(x,t) = —Ky{x,t)=^Jl^nn^du;diy (2.13) 

Another commonly used simplification to determine the proportion of extinction and 
emissivity by a source function is to impose local thermal equilibrium (LTE) where 
the radiation field is characterized by one single variable, the temperature T (black body 
radiation). In this case, the specific intensity depends only on the gas temperature and 
is given by Planck's law 

, 2hu^ 1 

/7,lte(x, t; u) = (2.14) 



e 



kT - 1 



and the collective, frequency integrated form yields the source function (x, t) : ]R^+^ 



^Many authors use a notation, where the frequency v is written as subset, instead of in the list of 
arguments /^(x, t; n) and the frequency integrated variables lack of that subset I,x, ■ ■ ■ ■ 
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5 1,4 



5, = — = y/,,LTEd.= ^T. (2.15) 



With X7,P we introduce a Planck mean opacity X7,P respectively the Planck mean opacity 
per mass k^^p 

oo 

/ X7-^7,LTE dv 

K^,PP = X-t,P = ^5 (2-16) 

/ -f7,LTE du 


and with X7,/? ^ Rosseland mean which will emerge in the radiation energy equation. 



1 _ i dT 

At this point we have defined all necessary variables and motivated a list of simplifications 
for our system of radiation hydrodynamics equations. However, two important steps are 
to be provided yet. First, one has to deal with the problem that in particular equations 
(2.9) but respectively all n moment equations contain in each case n + 1 momenta of 
the radiation field which is known as closure problem. The last step to the full set of 
equations of RHD will be the Lorentz transformation into the fluid frame. 



2.2.3 Radiation Closure Relation 

In section 1.1 it was mentioned, that imposing hyperbolicity on the conservation law can 
help to narrow possible solutions in some way. A conservation law is called hyperbolic if 
the Jacobian matrix associated with the fluxes Vjf has real eigenvalues and there exists 
a complete set of corresponding eigenvectors. As we have seen, interaction with matter 
adds source terms to our radiation conservation laws, we abbreviate them from now on 
with and S;^. 



dtJ-y + div (cF^) = J [r]^ — X'y^'y) du 



9tF^ + div (c^P^) = J [r]^ - X'yl'y)nduj (2.18) 



Si 



However, these terms do not influence if the system is hyperbolic or not since they appear 
without derivatives. We define the two Eddington factors (actually they are non-negative 
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tensors) 

-F=^, V = ^ (2.19) 

which can be interpreted as normahzed flux and pressure and therefore are a measure of 
the geometry of the radiative field. It is reasonable to assume that V is only a fuction of 
J-" and henceforth we write dV/dJ- = V' . Likewise we define the norms of the Eddington 
factors / = and p = \V ■ J^/f\- Now we express the density function and the fiux of 
the conservation law by means of the Eddington factors which represents our choice of 
dependend and independend variables for the Jacobian. 

= {Jj,cTJ^), f-y = {cTJ^,c^VJj) (2.20) 

The Jacobian of the radiative conservation law yields^ 

For upper Jacobian we obtain the following characteristic equation 

A = 2 (cp' ± V'^^ (4p + p/(_4/+p/))^ 

and system (2.20) is hyperbolic if the discriminant is not negative i.e. {p' — 2/)^ + 
4 (p — /^) > respectively as long as the constraint f'^<p holds. Moreover, / and p as 
normalized moments obey {p, /) < 1 and we derive the important conclusion 

f<P<l. (2.22) 

Equation (2.22) can be interpreted as flux limiter, which ensures causal values for the 
fiux function, namely that it can not get greater than the energy density J-y times the 
speed of light c. Physically it is clear that / and p take the value 1 for transparency 
which is often called the free-streaming limit. In this case, radiation propagates freely 
and the total amount of intensity is transported by the fiux function. The characteristic 
speeds in this regime evidently yield A = ±c. 

On the other end of the opacity scale, the optically thick regime, photons have very 
small mean free paths and underlie random-walk processes through the material. In 
this case, the radiation field is totally isotropic, the net flux vanishes = and the 
radiation stress tensor = gets diagonal (for arguments see Appendix 6.1.3). 

Consequently the eigenvalues of the Jacobian yield A = ±C\/l/3 which are the charac- 
teristic propagation speeds in this case. 



eusing ^ 



= and = -^rdF. 
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When we will discuss numerical high resolution methods for such nonlinear conserva- 
tion laws in section 3.1.4 we will encounter the concept of flux limiting again. To learn 
more about possible closure relations associated with the hyperbolicity approach, see 
Pons, Ibanez, Miralles (2000) [ ;5]. 



2.2.4 Fluid Frame 



In order to account for the local radiation properties of a fluid, we need to consider the 
interaction of matter and radiation in a common frame, ideally the comoving frame of 
the fluid flow. The transformation however not only has to describe advection effects 
(i.e. change in photon number density), but also impacts on momentum, energy and 
directions due to the Doppler shift of the photons. For the purpose of simplicity we will 
ignore pure scattering terms and assume local thermal equilibrium (LTE) again. 

An adequate procedure, especially with possibly relativistic fluid flows allowed is to 
apply Lorentz transformations. Introductions to Minkowski geometry and Lorentz trans- 
formation can be found in several appropriate text books like Ray d'lnverno [10]. At this 
point we just note that the relation between a four-vector x = (ci, r) in the lab frame 
and the moving frame x — >■ xq is given via 




7 -?u \ '^^ 

l + (7-l)^ ; ^ r 



2u 



(2.23) 



respectively x°o = L'^px^ with (a,/3) = 0,1,2,3. The subset denotes variables in 
the comoving frame, the gamma factor is defined by = yjl — |up/c^ where u is 
the fluid relative velocity. When we apply these rules to the photon four-momentum 
Pi = ^(l,n"), we obtain the required relations for the photon propagation in the fluid 
frame. 

/ n • u 

uq = 7Z^ ( 1 — 

no 



vq\ c 7-1-1 

These formulae for redshift (respectively blueshift) and aberration often are yet simplified 
for practical astrophysical applications, arguing that terms C'(^)^ can be neglected since 
fluid velocities are at least one magnitude smaller than the speed of light. This reduction 
also implies that the outstanding transformation is mere spatial as with 7 = = 1 
follows t = to- 

When we intend to transform the radiative quantities like the spectral intensity and its 
derived variables we need to occupy ourselves with relativistic kinetics respectively with 
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a relativistic form of the Boltzmann transport equation. To wit, the spectral intensity 
itself is not Lorentz invariant but the photon distribution function /y defined earlier. 
These considerations date back to Thomas (1930) [ !'^] respectively Lindquist (1966) 
[29], who drafted a general relativistic formulation. 

The connection between the spectral intensity, emissivity and absorption in the lab 
frame and in the comoving frame accordingly yield 

-^7 = (^) ^70' ^7 = (^) ^70' ^7 = ^^70- (2-25) 

The frequency integrated moments of the radiation field occur as elements of the radia- 
tion energy-stress tensor which is defined by 

Tf = c' f %y^f, (2.26) 
and can easily be rewritten to a more familiar form''. 

The conservation law associated with radiation transfer gains its most elegant form in 
this notation, namely V • T with source terms S-y = (5^,8;^) as defined before. 

The Lorentz transformation of yields T-y^ = LT-yL as the operator is unitary and 
thus L = L^.s When we confine ourselves to small velocities again and discard terms 
higher 0(^)^, we would receive the applicable approximation formulae. 

What is left to do, is to mark the moments of the transport equation in this classical 
limit. For forcible deduction we refer to Buchler (1979) [ ] and (1983) [ ] who elaborated 
a fluid frame description in arbitrary geometry inclusive of terms C'(^). The frequency 
integrated radiation energy equation and the radiation flux equation in the comoving 
frame accordingly yield 

dtJyO + div (u Jyo) + c div H.yo + K^o : grad u - cXj,po{Jjo - S^o) = 
(9tHyO + div (uH^q) + c div K^o + H^o • grad u + cx7,ijoH^o = 0. 

(2.28) 

Notation: from now on, the neat radiation variables ( J, H, . . . ) shall be read as photonic, 
frequency integrated functions in the comoving frame {J^QjU^o, . . .). 



^Using = ^(l,n") and = hu. 
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2.3 Gravitation 

In section 1.2 we highlighted some mathematical basics for the Laplace and the Poisson 
equation in context of elliptical partial differential equations. The Poisson equation for 
the gravitational potential, a scalar function <I>gi.av(x, t) : M"^^ — ^ M yields 



where /?ni(x,t) : W^~^^ — ^ M"*" is the time dependend mass density again in n = 1,2,3 
spatial dimensions. 

However, most applications of radiation hydrodynamics with gravitation avoid solving 
the Poisson equation for the gravitational potential per se, especially in spherical sym- 
metric applications'^. The main challenge associated with self gravitation in 2D and 3D 
is that it requires the numerical solution of the Poisson equation in each time step which 
is computationally expensive. Moreover, it is not trivial to find appropriate boundary 
conditions for this elliptical problem in astrophysical context since most objects are not 
bounded naturally but artificially respectively numerically. 

2.3.1 Self-Gravitation in 2D and 3D 

There have been efforts in finding and applying solution methods for self-gravitating 
flows in multiple dimensions since decades. Differential algorithms solve the Poisson 
equation in the form (2.29) whereas integral methods rather consider the theory of 
Greens functions as briefly mentioned in section 1.2 where the formal solution in 3D 
yields 



With latter family of algorithms the grids boundary causes no problems and the inte- 
gration could principally be conducted with standard methods. However, the number of 
grid cells to be summed (~ A^^) and therewith the computational cost get huge in mul- 
tiple dimensions. Several authors thus suggest to expand the integrand e.g. in spherical 
harmonics in order to reduce operations to an order ~ N. Miiller and Steinmetz (1995) 
[34] developed an efficient Poisson solver based on such a separation ansatz. Cosmo- 
logical hydrodynamic simulations have affiliated these ideas and implemented them in a 
number of cosmological codes. 

Spingel (2009) [39] suggests an interesting method to realize self gravitation in a con- 
servative way by reformulating the total energy equation. He adds a gravitational source 

*With radial symmetry, the gravitational acceleration G(r) is simply determined by the integrated 
mass m(r) = ^'Kpm{r')dr' via the well known relation G(r-) = — Sr'i&grav = — , an integrated 
form of the Poisson equation. 



A<I> 



grav — 



'm 
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term on the right hand side of equation (2.5) whereby it gains the form of a conservation 
law that ensures the total energy to remain constant. The total energy of the system 
consequently yields the sum over all integrated cell energies. In stars however this mod- 
ification is rather inapplicable since the contribution by potential energy exceeds all 
other energy terms by far and this term would numerically dominate the energy equa- 
tion which is unfit in practice. Hence we stick with the formulation where changes in 
potential energy directly contribute to the momentum equation (see the full set of equa- 
tions 2.4). Alternatively we suggest a non-static equation for the gravitational force in 
multiple dimensions. 



2.3.2 Non-Static Gravitation Equation 

The gravitational potential <I>grav is not needed explicitly at any point of the computation 
since coupling is accomplished via the gravitational force respectively the gradient of the 
potential i.e. grad^grav 

We define the gravitational force vector G(x, t) : M"^"'^ — ^ M"+-'^ in the familiar way 

G = -grad$grav (2.30) 

hence the Poisson equation can be written div G = —AirGp. We differentiate with respect 
to t and use the continuity equation (2.2) to obtain 

dtdiwG = -A^Gdtp 

= AttG div (pu) . 

Since $grav is C°° it is reasonable to assume that the derivatives with respect to t ex- 
ists and is continuous as well, the operations can be interchanged. We consider the 
conservative, integral form of these equations again and use the divergence theorem. 



j diwdtGdV = J AirG div (pu) dV 

V{t) Vit) 

J dtG-dS = AttG J {pu) ■ dS (2.31) 

dVit) dV{t) 
We obtain an equation for the gravitational force in form of a pure initial value problem. 

dtG = 47rGpu (2.32) 

This way we avoid solving the Poisson equation in multiple dimensions by solely con- 
sidering its gradient, the gravitational force. Equation (2.32) describes gravitation non- 
statically by understanding temporal changes of gravitational forces as mass transport. 
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If, respectively in what extent, this formulation is numerically preferable will have to be 
studied in future computations. 

Self-gravitation expressed via (2.32) would not have to be solved at any time of the 
evolution but only initially contrary to Poisson equation implementations. Hence it is 
supposedly especially suited for explicit and parallelized methods. 
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2.4 Full Set of RHD Equations with Self-Gravitation 

The system of equations of radiation hydrodynamics with gravitation is now presented 
in its for our purposes somewhat simphfied formulation. The fohowing system has been 
basis for a number of imphcit RHD computations hke Dorfi (1999) [ ] and Stokl (2006) 
[40]. In latter PhD thesis physical assumptions and simplifications are illustrated and 
their consequences executed more detailed than in this paper. 

Continuity Equation 

dtp + div {up) = (2.33) 

Equation of motion 

dt{pu) + div (puu + P-a) - pG- ^Xi?H = (2.34) 

Equation of Internal Energy 

dt{pe) + div (peu + (P - o") • u) - 4ttxp{J - S) = (2.35) 
Equation of Radiation Energy 

dtJ + div (uJ) + cdivH + K : gradu - cxp{J - 5) = (2.36) 
Radiation Flux Equation 

dtU + div (uH) + c div K + H • grad u + cxrR = (2.37) 

Poisson Equation 

= A-kGp (2.38) 

respectively 
Non-Static Gravitation Equation 

5tG = 47rGpu (2.39) 

Basically we collected equations (2.2), (2.4), (2.5) plus (2.28) and added gravitation. 
The viscosity a will be replaced by an artificial viscous pressure — Q which we establish 
in section 3.1.4. In section 3.3 we present this system of equations in strong conservation 
form for non-steady curvilinear coordinates. 
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In this chapter we want to develop some numerical methods for conservation laws after 
having motivated and anticipated relevant mathematical concepts in the previous sec- 
tions. Concepts for error functions, convergence and stability for linear and nonlinear 
equations will be studied in form of an overview. A particular focus will lie on conserva- 
tive methods for nonlinear problems in general curvilinear coordinates which will require 
some tensor analysis and differential geometry. 

In sections 3.2 and 3.3 we will present the mathematical framework for strong conser- 
vative numerics in non-steady curvilinear coordinates and exemplarily sketch the system 
of RHD in spherical and polar coordinates. An articulate introduction and motivation 
to computational methods in astrophysical context can be found in LeVeque, Mihalas, 
Dorfi, Miiller [25]. 

3.1 Numerical Methods for Conservation Laws 

Standard numerical methods for partial differential equations are established under the 
assumption of (classical) differentiability. Routine finite difference schemes of first order 
usually smear or smoother! the solution in the vicinity of discontinuities, since they come 
with intrinsic numerical viscosity. Standard second order methods show something to 
the effect of the Gibbs phenomenon, where we see oscillations around shocks. In the 
past decades so called high resolution methods have been developed in order to achieve 
proper accuracy and resolution for nonlinear, discontinuous problems as they appear 
also in radiation hydrodynamics. In the following sections we will introduce some basic 
technical terms in discrete numerical methods for conservation laws and present some 
important concepts for nonlinear problems. 

3.1.1 Introduction to Discretization Methods 

The basic idea of discretization methods for PDEs is to represent functions by a finite 
amount of data and to approximate differential operators by difference operators. To 
simplify matters, we confine ourselves to 1 -|- 1 dimensional problems and constant in- 
crements in the upcoming introductory sections. Our toy model for analysis shall be 
the time-dependent Cauchy problem (1.5), namely 9td -|- div f (d) = but in linear form 
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and one space dimension^ dtd + A dxd = with constant A. Concerning structure and 
nomenclature we basically follow LeVeque [20] in the following sections. 

For the spatial part we introduce a mesh width h = Ax and a time step k = At for 
the temporal discretization. Any point in this computational domain can be attained 
by multiple steps in t and x, that is {xj,tn) = {jh,nk), where j G Z and n G N with 
starting point or initial condition d{xj,to)- What we want to obtain is the approximate 
solution of the exact d{x,t) at discrete grid points d{xj,tn) and we characterize these 
approximation as D". 

The difference scheme is somewhat arbitrary and its design depends on the concrete 
problem to be approximated. The essential requirement is that in the limit of (j, n) — >■ 
the difference operator reproduces the differential operator. The method is called explicit 
when the function at the new time step is calculated directly via function values at prior 
time steps. We introduce a special notation for such operators with T-L(D^;j) = D^"*"^ 
which reproduces the vector of approximations at the new time step D"^^. Some well 
known explicit schemes for our toy model are for example the so called Upwind, Lax- 
Priedrichs, Leapfrog and Lax-Wendroff method. 



k 



?^(D";i) = B]--A{B]-B]_ 



-1) 



1 h 

2h 

k 



-(D^_,+D,V)-^^(Di+i-Di-i) 



(3.1) 



3.1.2 Error and Stability in the Linear Case 

There is a number of criteria that reveal the quality of a numerical method, however the 
most obvious is to define an error function which states the difference between the exact 
and the numerical approximate solution. In this sense the pointwise error is defined by 

E" = D" - d" (3.2) 



^We maintain the tensor notation in boldface d for densities and fiuxes for distinction purposes als in 
1+1 dimensions. 
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In conjunction with conservation laws the numerical solution is often seen as an approx- 
imation to a cell average of d(xj, tn), defined by 



d] = ^ f d{x',tn)dx', (3.3) 



i-1/2 

where the half step is simply Xj_|_i/2 = ~\~ h/2 and the pointwise error in terms of cell 
averages yields 

%" = - d^. (3.4) 

Let be Dfc(x,t) = D" for (x,t) G [a:^j-i/2; ^j+1/2) ^ [^min+i), then the error function in 
terms of piecewise constant functions is defined by 

Efe(x,t) = Ufc(x,t) -u(x,t). (3.5) 

We have already mentioned that the difference operator has to approximate the differ- 
ential operator properly. Regarding the temporal iteration we want that smaller time 
steps cause smaller errors. We call the numerical method convergent (in a norm'^), if for 
arbitrary initial condition do, Vt > 

lim||Efc(.,t)|| =0. (3.6) 

fc— 5-0 

Assuming smooth solution we always can expand the exact solution into a Taylor series 
and compare the leading terms to the numerical scheme. By that means we get a measure 
of how well the difference equation reproduces the difference equation locally. The local 
truncation error is therefore defined by 

U{x,t) = ^(d{x,t + k)-U{d{.,t)-x)) (3.7) 

and we call the numerical method consistent, if 

lim||Lfc(.,t)|| =0. (3.8) 
fc— >o 

The method is of order p, if for all sufficiently smooth initial data with compact support, 
there is some constant Cl, such that 

||Lfc(.,t)|| <ClF, (3.9) 



^For conservation laws the a practicable norm is the 1-norm || . . . || = 1 • • • 1 dx, since more rigorous 
demands come into conflict with approximation of discontinuities; ||Dfe(., t„) ||i = . 
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for all k < ko and t < T. The restriction to a local T falls for stable methods, where the 
order remains the same for any t ox k. 

There is a profound stability theory for linear methods; its details can be found in 
several relevant text books (e.g. Richtmeyer, Morton [•')()]). Nevertheless, at this point 
we shall mention the Lax- Richtmeyer stability which states that a method is stable if for 
each time T there is a constant Cs and a value /cq > such that 

\m\\ < Cs (3.10) 

for all nk < T, where k < ko. The important Lax equivalence theorem claims that for 
consistent linear methods stability is necessary and sufficient for convergence. 

In this context we also should mention the Courant, Friedrichs, Lewy ( CFL ) condi- 
tion, which seizes the fact that hyperbolic equations stand out due to a bounded domain 
of dependence. The solution at a point d(x, t) is determined by initial data within some 
finite distance to this point. This range of influence is related to the fact, that hyperbolic 
equations always have finite propagation speed, as briefly mentioned in the introduction 
to section 1.1. When Courant, Friedrichs and Lewy studied convergence of numerical 
schemes for hyperbolic conservation laws, they observed that one necessary condition for 
convergence is that the domain of dependence of the finite difference scheme includes the 
domain of dependence of the PDE. In other words, data from one time point must not 
propagate faster than the characteristic speed of the original equation. Mathematically 
worded, the domain of dependence V{x,t) is a set of points {x,t) with the property 
that a disturbance of d(x,t) influences the solution at (x, t). The numerical domain of 
influence Pfc(x, t) is analogously defined. For the satbility of an (explicit) method it is 
necessary that for /i, A: — >■ the numerical domain of dependence ^ifc(x, t) contains asymp- 
totically the domain of dependence the corresponding initial value problem. Practically, 
an inevitable (but not sufficient) stability condition for such methods is represented by a 
limited ratio of time step to mesh size. For a scalar advection equation with propagation 
speed a, the CFL condition yields > ^ > 1 as an example. 

3.1.3 Conservative Methods and Lax Wendroff Theorem 

In the linear case, the Lax equivalence theorem states that a consistent and stable dis- 
cretization scheme automatically converges to the weak solution of the conservation law. 
For nonlinear equations, consistent methods do not necessarily converge to the weak 
solution of the conservation law'^. We have already mentioned that nonlinear conserva- 
tion laws can have several weak solutions but we are only interested in those, that also 

^E.g. the CIR method, a nonlinear generahzation of the Upwind scheme, is not applicable for discon- 
tinuous solutions, since it is not conservative. 
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satisfy the entropy condition. So it is not surprising that it is required to make sure, the 
numerical scheme satisfies some discrete analogon of the entropy condition (1.14). 

The proper postulation to the numerical scheme that prevents convergence to non- 
solutions is that the difference method has conservation form, i.e. 

= D« - ^(f(D^_,, . . . ,D^^,) - F{U]_^_„ . . . ,D^^,_0) (3.11) 

for a numerical flux function F of p + g + 1 arguments. This numerical flux can 
be interpreted as averaged flux through Xj_^_i/2 and [tn,tn+i] in the sense that F = 
^ /j*^'"''^ f (d(xj_|_x/2) 0) In section 3.2 we will turn our intensive attention to this 
numerical quantity when we concern ourselves with conservative methods in general 
coordinates. 

We also reformulate the concept of consistency for conservative schemes. The con- 
servative method is consistent with the original conservation law, if the numerical flux 
function F reproduces the true flux f in the case of constant flow d so that 

F(d,...,d)=/(d). (3.12) 

Moreover, we require as the arguments of F approach a common value d, that the value 
of F meets f(d) smoothly^. So for each d there may exist some constant > so that 

F(Dj„p, . . . , Bj+q) - f (d) < K max Bj+i - d . (3.13) 



Before we can establish a nonlinear convergence theorem we first have to define a total 
variation of an arbitrary function v. 

oo 

TV(v) =limsup- / ||v(x) - v(x-e)|| (3.14) 



-oo 



which reduces to a more familiar look for differentiable v respectively appropriate dis- 
tribution derivatives, namely TV(v) = ||v'|| dx. We call a sequence of solutions D, 
convergent to a solution d, if for every bounded set [a,b] x [0, T] 

T b 

\\Di{x,t) -d{x,t)\\ dxdt for i ^ oo (3.15) 

0. 



*The requirement is Lipschitz continuousity which is exigently given for any smooth and consequently 
differentiable flux function f . 
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and there exists an upper bound R for the total variation such that 

TV(D,(-,t)) < (3.16) 

for < t < T and i G N. 

This leads us to the essential Lax- Wendroff theorem which averts the following. Let 
{Di(x,t)} be a sequence of solutions of a conservative and consistent difference method 
for the one dimensional conservation law, where the mesh parameters {k, /i) — ^ for 
i —7- oo. Suppose the numerical approximation at the i-th grid {Di{x,t)} converges to a 
function d(x, t) in the sense made precise before; then d(x, t) is a weak solution of the 
conservation law. 

Still, the consistent conservative scheme can converge to a weak solution that does not 
satisfy the entropy condition (1.14). We add a comparable condition for the numerical 
method in form of the Lax- Wendroff corollary. Let a entropy function and (p entropy 
flux and let be $ the consistent (in the above defined sense) corresponding numerical 
flux function. Given the Lax- Wendroff theorem is satisfied and 

a(D]+') < a(D]) - ^(^(D",i) - cI>(D", j - l)) (3.17) 

holds, then d(x,t) satisfies entropy condition (1.19). 

As in the linear case, the picture is not drawn completely until we clarify stability of the 
numerical method. All we know now is, that if a sequence of approximation converges, 
the limit is a weak solution. However, we require intrinsic stability of the method in 
order to ensure convergence of the difference scheme. Without going into details, total 
variation stability is the suitable concept for nonlinear methods, for which the Lax 
equivalence theorem is not applicable. It can be shown that TV-stability, consistency 
and conservativeness guarantee nonoscillatory convergence of the advection scheme. 

As a further remark to the topic of methods for nonlinear conservation laws, we want 
to present an important class of methods that ensure TV-stability. The idea goes back 
to Harten (1983) [IS], who was studying high resolution methods for hyperbolic conser- 
vation laws and proposed so called total variation diminishing (TVD) (originally non- 
increasing) methods. We call a numerical scheme 7i(D'^;j) total variation diminishing 
if 

TV(D"+i) < TV(D") (3.18) 

for all grid functions D*^ at all time steps n. Of course, the true solution of the conser- 
vation law must be compatible with this postulation and indeed it can be shown, that 
any weak solution d{x,t) satisfies TV(d(-, t„+i)) < TV(d(-,tn)) for any n > 0. 

In the upcoming two subsections we want to present some sophisticated high resolution 
numerical methods incorporating these techniques for nonlinear conservation laws. The 



36 



3.1 Numerical Methods for Conservation Laws 

emphasis will lie on concepts and techniques relevant for our RHD problem as they have 
been implemented in applications like Dorfi (1999) [ ] respectively Dorfi et al. (2004) 
[12] 

3.1.4 Artificial Viscosity 

The TVD approach is not the only one to find entropy-satisfying weak solutions with 
numerical schemes. One of the simplest techniques are monotone methods"", however the 
have one main disadvantage. Monotone methods are at most first order accurate. In 
the past decades major efforts were made in developing numerical methods for nonlinear 
hyperbolic conservation laws that are at least second order. 

One patent attempt to finding such a high resolution method is to adapt a well known 
high order method like Lax-Wendroff for nonlinear problems. As we have contemplated 
in section 1.1.3 we can add an artificial viscosity term to the conservation law in a way 
that the entropy condition is satisfied. Numerically we are keen to design this viscosity 
in a manner that it affects discontinuities but vanishes sufficiently elsewhere so that 
the order of accuracy can be maintained those regimes where the solution is smooth. 
Interestingly this idea was inspired by physical dissipation mechanisms and dates back 
more than half a century to VonNeumann and Richtmeyer (1950) [47] who investigated 
new methods to treat shocks in hydrodynamic simulations. Bearing in mind that we 
have to provide conservation form, the numerical flux function gets modifled by a an 
artiflcial viscosity q(D; j) for instance in the following way. 

Fvi.c(D; j) = F(D;i) - /iq(D; j)(D,+i - D,) (3.19) 

Since the original design of this additional viscous pressure in the form q = C2/5(Au)^, 
C2 G M, for one dimensional advection dtd + adxd = qdxxd it underwent a number of 
modifications and generalizations. It turned out to be numerically prerferable to add a 
linear term, see Landshoff (1955) [ ] in order to control oscillations, generalizations to 
multi-dimensional flows mostly retained the original analogy to physical dissipation and 
reformulated the velocity term accordingly, e.g. Wilkins (1980) [49]. 

In any case artificial viscosity broadens shocks to steep gradients at some characteristic 
length scale but on the other hand shall not cause unintentional smearing. The concrete 
composition and implementation of this artiflcial viscosity coefficient q depends on the 
application. Tscharnuter and Winkler (1979) ['■ '■] pointed out that the viscous pressure 
in 3D radiation hydrodynamics has to unravel normal stress, quantified by the divergence 
of the velocity field and shear stress, which is expressed by the symmetrized gradient of 
the velocity field according to the general theory of viscosity. It is designed to switch on 

^The more general concept is called Zi-contracting methods, see e.g. Harten et al. (1976) [19]. 
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only in case of compression (div u < 0) yielding 



Q = — g|/^igj./3max(— div u, 0) ( Vu — -edivu 



(3.20) 



and is set zero otherwise. In section 3.2.4 we will elaborate, why we have to modify this 
definition for general geometries by replacing the unit tensor e by a metric tensor g. In 
section 3.3 we will specify co- and contravariant components of the artificial viscosity 
tensor in spherical coordinates. 



3.1.5 Flux Limiters 

When we discussed the radiation closure relation in section 2.2.3, the approach was ana- 
lytical respectively physical. Hyperbolicity of a conservation law implies real eigenvalues 
of the corresponding Jacobian Vdf which leads to characteristic propagation speeds 
which in case of radiation are connected with the speed of light. Pons, Ibanez, Mi- 
ralles (2000) [35] studied flux limiting in context of Gudonov-methods based on these 
constraints on propagation speeds, Levermore (1984) [27] studied flux limiters and Ed- 
dington factors in anisotropic cases. 

Flux limiters (or slope limiters) basically restrain solution gradients near shocks and 
thereby avoid oscillation in high resolutions methods (such as second order TVD). With 
flux limiting, the advection scheme gets split up in high order fluxes in smooth regions 
and lower order schemes (like Upwind) in the vicinity of discontinuities. With a limiter 
function $(D; j) the flux yields F = Fiow + ^{F high ~ Flow) respectively in a more 
pracitical formulation 



There is a number of different ways in implementing such a flux limiter albeit they are 
typically functions of the steepness of the gradient. One ambiguity lies in this measure of 
smoothness in discrete terms, in this context commonly denominated 9j. One passable 
way is the ratio of consecutive gradients 



and is obviously near 1 for smooth data. Several flux limiter methods, i.e. different 
functions ^{6j) have been suggested and tested extensively (e.g. Van Leer (1979) [15]). 
In implicit RHD computations like Dorfi et al. (2004) [12] and previous papers implement 



F(D;j) = Fhigh(D;i) - (1 - $(D; j))(Fhigh(D; j) - Fiow(D;i)). 




9 



D,+i-D, 



(3.22) 
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0.0 0.5 1.0 1.5 2.0 2.5 3.0 

Figure 3.1: Van Leer Limiter 



the Van Leer limiter which is constructed in the following way. 

29 



1 + 




if 61 > 
otherwise 



(3.23) 



The Van Leer hmiter describes a smooth function, a somewhat distinct characteristic in 
comparison to other feasible methods. In figure 3.1 the gray region denotes TVD; for 
details we refer to LeVeque [_ ] once again. 



3.2 Conservative Methods in General Curvilinear Coordinates 

In the previous section we presented necessary premises for conservative difference schemes, 
particularly the conservation form for a numerical method (3.11). The numerical flux 
function F we defined can be interpreted as flux averaged over a cell (3.3) and a time 
interval F = ^ It^^^ f (d(^j+i/2) ^)) dt. Descriptively written, the conservation form en- 
sures that whatever fraction of the density gets transported from one to another cell, 
arrives there. In other words, if supposedly all density of one cell gets advected to an 
other, we need to find the same physical amount in there afterwards. However, if the 
cell volumes differ and the flux is a cell integrated quantity it is required to assess the 
cell volume itself properly. Hence in order to conserve a certain density function d(x, t) 
in general coordinates where non non-steady, nonlinear volume elements occur, we need 
to investigate the numerical flux more thoroughly in terms of differential and integral 
calculus. 

Mathematically the fluxes that occur in hyperbolic conservation laws are generated via 
spatial differential operators which are obtained by applying the divergence theorem to 
a differential volume, bounded by a surface (1.4). In the upcoming subsections we want 
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to elaborate how differential geometric considerations lead to a reformulation of these 
differential operators, Gradient and Divergence in a geometrically conservative form. 



3.2.1 Tensor Calculus and Elementary Differential Geometry 

To begin with, we want to recapitulate some relations from tensor analysis and elemen- 
tary differential geometry and define linear and nonlinear coordinate systems, coordi- 
nate transformations, base vectors and co- and contravariant vector components. In the 
following section we will use Einstein summation convention, indices in braces do not 
designate components but coordinate sytems. Consequent argumentations and proofs 
can be found in Lichnerowicz [28], Schutz [-IT] and Liseikin [.M]. 

Let A4 be an affine space with corresponding vector space V and S(q,) a linear coordi- 
nate system. The coordinate system maps points x £ M (linear, bijective) to a d-tuple 
of numbers x*. : Ai — t- M'^, x ^ (x^q,)*). Coordinate lines of a linear coordinate 

system are straight lines on the affine space. Two linear coordinate systems '^[a) s-^d 
S^^-j are connected via a coordinate transformation 

X(a)' = A(„^)j X(^pf + A(„^)\ (3.24) 

Such affine transformations imply translation (A(fj^)*) as well as rotation and shearing 
(A(q,^)P of coordinate lines. The natural base vectors of any coordinate system are 
tangential vectors to the coordinate lines and thus constant in the linear case. 



^ (3.25) 



The index i in this case denotes the i-th base vector. The dual base is implicitly defined 
via the inner product ©(q)* '^{a)j = ^^j] t^ie dual base e(^)* G V* is a linear map V* — )■ M, 
which maps the base vectors e^^-jj £ V to zero or one. The natural and the dual base 
vectors transform unequally to a new coordinate system Sj-^). 

^W=^o.,)]e^»)' = Q^^io^)' (3-26) 

The matrix ^(ai3) is referred to as Jacohian of the coordinate transformation. 

A transformation to a curvilinear coordinate system is no longer necessarily linear 
nor bijective. With nonlinear coordinate systems, the base vectors are not constant but 
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functions of position. 

Of course, the entries of the Jacobian A/^g) are now functions too 



A(.^);.(x) = (3.28) 

and e(^)*(x) • e(Q,)j(x) = is vahd pointwise. All vectorial and tensorial quantities T;; 
in this domain are defined via their components with respect to the local base in the 
following sense. 

^Wp.K^) = ^(^W'W' • • • ,e(„)-''(x),e(„),(x), . . .e(,),(x)) (3.29) 

We define the covariant components of the metric tensor g as local measure of distances 
in the coordinate system 

ei(x) • ej(x) = ^^^(x). (3.30) 

The metric tensor is obviously symmetrical due to commutativity of the scalar product 
and implicates a symmetrical bilinear form x • y = gijX^y^ on the Euclidean vector 
space V. From now on lower indices designate covariant and upper indices contravariant 
components. Together with the contravariant metric tensor e*(x) • e-'(x) = g^^(x) we 
gain the following transformation rules between co- and contravariant components of a 
tensor. 

Vi = QijV^, v' = g'^Vj (3.31) 

Raising and lowering indices in linear coordinates just means to transpose respectively 
to interchange rows and columns of higher rank tensors. With curvilinear coordinate 
systems the dual space is also metrically different. 

Next we want to liaise the recently defined geometrical items with differential calculus. 
Vectorial and tensorial quantities feature a native construction in differential geometric 
terms. Vector components that are generated via a Gradient of a scalar are covariant by 
definition, e.g. di(f) = pi whereas coordinates are contravariant. Higher rank tensors 
are constructed by differentiation or tensorial products, e.g. t = pq*^, and its native 
depiction in components is ti^ = PiQ-' ■ The explicit form of such a tensor t in coordinates 



x' is 



t(x) = fi^(x)eXx)ej(x). (3.32) 



For this reason, the tensorial physical quantities in the equations of radiation hydrody- 
namics will have to be considered firstly in their native form as well. When performing 



^We generally omit explicit tensor products pq = p ® q respectively e*e-' = e* ® e-' . 
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operations at and between them it will be necessary to consider a consistent denotation 
(e.g. contravariant components) and raising and lowering indices according to (3.31) will 
be required. 

Evidently when differentiating a tensor we have to consider not only its components 
as functions but also the generally non constant base vectors in the coordinate system 
S(q,). Without forceful motivation we define a covariant derivative of a vector field '^(a)^ 
and the Christoffel symbols T(a) the following way. Let V be a map that transforms 
a tensor of rank (m, n) to a tensor of rank (m, n + 1), satisfying ^(p = d(j) for any scalar 
(j). Let it furthermore satisfy all congenital demands on a derivative like the Leibnitz 
rule, additivity, commutativity and torsion- freeness and it shall be associated with the 
partial derivative via 

Viv^ = div^ + T^kV^ (3.33) 
where the Christoffel symbols compensate for the locality of the base vectors. 

r^',fc(x)ej(x) = diBk{-K) (3.34) 

We postulate that this covariant derivative provides the proper afRne connection am im- 
plicitly defines the correct transformations between general coordinate systems and can 
be expanded to tensors of higher rank and arbitrary composition of co- and contravariant 
components. 

diV^ + T\kv'^ = div V 
diVj - T^ijVk 

ditj -\- r litj r jitl 

(3.35) 

For linear coordinate systems the Christoffel symbols vanish trivially and the covariant 
derivative is identical to the partial derivative. Moreover, the covariant derivative can be 
called metrical in terms of compatibility with the metric tensor g. The Ricci Lemma for 
the fundamental tensor states that the covariant derivative of the covariant components 
of g vanish, which can be checked quickly. 

Vi5jfc = digjk - T'-ijQik - T^kQji = (3.36) 

When we contract upper equation with g^^ we obtain g^^digjk — T^ij — T'^ik = and gain 
an interesting formula 

r^ifc = l9''d,gjk (3.37) 
which reminds of the derivative of the determinant of the metric tensor det g = g which 



= 
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IS 



di9 = gg^^diQjk- 



(3.38) 



This yields to the following important relation between the Christoffel symbols and the 
metric. 

1 dig _ Oi/|g| 

2 5 " 



J- ki 



(3.39) 



We recognize the left hand side as a geometrical term in the divergence of a vector field 
(3.35) which leads to a meaningful reformulation of this derivative. 



divv 



1 



-A 



|g| 



(3.40) 



Immediately we see the connection to the divergence theorem from integral calculus 



when we remember that the Jacobian determinant defines the volume element dV( 



(a) 



\^{a)\dx'^ and we know that |A(-cf)| 
div V dV = 



5(a) I 



-A 



|g| 



|g|^^ 



|g| dx^ . . . dx^' 



d. 



dx^ . . . dx" 



(3.41) 



Latter expression equals the flux of the vector field v over the boarders of the volume, a 
factor which shall be crucial for our geometrically conservative formulation of the equa- 
tions of RHD. Unfortunately the clue (3.41) is conflned to scalar equations which is 
expressed in Vinokur's theorem published by Vinokur (1974) [ IG] and proven elegantly 
by Bridges (2008) [5]. It states, how strong conservation form (i.e. geometrically con- 
servative in our diction) can be maintained regardless of what coordinate system used 
also for vectorial conservation laws. 



3.2.2 Strong Conservation Form 

The numerical flux function we introduced in (3.11) as integral quantity has to be un- 
derstood not solely as a function at nodes but as a function of cells and their geometries, 
that is to say 

F = F(D,^,e*). (3.42) 

Before we postulate how we obtain geometrically conservative fluxes or in other word, 
how we need to rephrase the covariant derivative (3.35) generally, we gain insight to 
another eclectic concept in differential geometry. The theory of differential forms is an 
approach to differential calculus in more dimensions that is independent of coordinates. 
Again we have to refer to relevant literature like Bishop and Goldberg [ ] for stringent 
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motivation and definitions. For us, differential forms shall just be objects that emerge 
under integrals, like wdx, fdS, (pdV, . . . and in that manner we have seen them consis- 
tently. We are going to like them because our (numerical) fluxes are objects that emerge 
under integrals as well. Moreover, we accept that there is a peculiar generalized con- 
cept of differentiating such differential forms which basically simplifies to the covariant 
derivatives when we settle for particular coordinates. 

We introduce a new notation, the wedge product of two vectors p and q as the following 



antisymmetric tensor^ 



pAq = pq — qp. (3.43) 



Tensor analytically this wedge product is motivated from the observation that any an- 
tisymmetric tensor t can be decomposed to linear combinations t = f^{eiej — ejei), 
{i < j) respectively in terms of the wedge product t = f^ei A e^. 

We define a p-form (i.e. a differential form of degree p) by 

w = Wi^,„ip dx'^ A dx^^ A ... A dx^" (3.44) 

which is obviously totally antisymmetric in its indices. The exterior derivative doj of a 
j»-form u) = Wi-^„,ip dx^^ A ... A dx^p is 

duj = djWi^,„ip dx^ A dx"^ A ... A dx"" (3.45) 

which yields a (p-l-l)-form. As an example, the exterior derivative of a 1-form uj = Wi dx^ 
yields duj = diWj dx^dx^ . The general Stokes theorem shows the connection of the exterior 
derivative on differential forms to integration theory, namely 



j u. (3.46) 



n an 

The generalized concept of the covariant derivative maps a p-form to a (p — l)-form and 
is defined by 

6uj = *d*uj (3.47) 

where * is the Hodge star operator which could be understood as generalization of the 
Levi-Civita tensor e with metrical weight. 

* (dx^i A ... A dx'p) = dx'p+^ A dx'p+^ A ... A dx'" (3.48) 

The Hodge operator is the source for the volume form wvoi = *1 = •\/[g|dx-'^ A ... A dx" 
which we already encountered in the divergence theorem (3.41). 

The benefit of dealing with differential forms shows when we associate them with well 



^Of course we imply tensor products pq = p ® q again. 
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known concepts of tensor calculus. The exterior derivative of an 1-form uj in yields 
the Curl, to wit d{'w dx) = curl w dS which in component notation is often expressed by 
d{w dx) = diW^ EijkdS^ . Moreover, as pointer to the nativeness of dealing with differential 
forms we mention that we mathematically defined the Curl of a vector field via a closed 
line integral (in a plane orthogonal to a normal vector n) in the limit of vanishing path, 
i.e. (curlcbi) • n = lim^_!.o ^§ oj • dr^. The main difference to tensorial concepts is, that 
with operations on differential forms we map something that emerges under an integral 
to something that emerges under an integral. Without further excursions we claim that 
we now have the munition to design geometrically conservative differential operators for 
general coordinates. 

The first elaborate insight to geometrically conservative differential operators can be 
found in Thompson, Warsi, Mastin [ ], although they completely exclude differential 
forms and play extensively with tensor analysis what makes the book a little delicate 
to comprehend. However, differential form theory perfectly endorses their conclusions 
and also Bridges' proof (2008) [ ] of Vinokur's theorem resorts to differential forms as 
virtually inherent description for conservation laws. We conclude the connection between 
covariant derivatives in terms of (3.35) and differential forms. The covariant derivative 
on a 1-form u = Wi dx^ yields explicitly 



which we recognize as result (3.40), namely divw. For the strong conservation form 
of the equations of radiation hydrodynamics we will basically need the following three 
operations. 

Gradient of a scalar </> 



'Let he Lj = a dx + b dy + c dz, then dw = {dyady + dza dz) f\ dx + ■ ■ ■ = {dyC — dzb) dy A dz + {dza — 
dxc) dz A dx + {dxb — dyo) dx A dy = rota; where dy A dz = , dz A dx = etc. 



5uj 



*d* LO 




(3.49) 




(3.50) 
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Divergence of a vector field v 



div V = V,?;* 



-A 



e • V 



(3.51) 



Divergence of a second rank tensor t 



divt = Vif^ej 



e'-t 



(3.52) 



Of course these differential operators can be translated to the well known covariant 
derivatives (3.35) with tensor analysis as well. 

1 



|g| 



r^^e^0 + a^(5^%)(/. + e^a^^ 



(3.53) 



The canonic shape of tensorial conservation laws / (5jd + divf(d)) dVdt = in strong 
conservation form for non-steady coordinates'^ yields 



dt 



|g| d 



f (d) • 







(3.54) 



where of course d and f still have to be decomposed according to their native tensorial 
components. 

We substantiate upper proposition by an example of a vector conservation law, the 
equation of motion (2.4) in n dimensions. The differential expression 



dt 



pu +di 



{puu + P + cr) • = 0, 



1, 



, n 



(3.55) 



implicates the integral form of conservation of momentum which would be considered 
numerically for non-steady coordinates. The important difference to a component-wise 
structure is that in this case undifferentiated terms would arise in n equations (see 
Vinokur (1974) [ ]). 

Neither the covariant derivative form for the divergence V(.) = d{.) +r(.), nor the k- 
th component of the expression 5j[-\/|g|/*'^] which would even be analytically true only 
for antisymmetric tensors f anyway, satisfies the conservation law adequately. Since 
the base vectors e* are functions of position for general coordinates, a separation as 
implied would introduce undifferentiated terms. However, the k Cartesian components 



^In the following section we will argue, why the volume element ^/\g\ also appears in the time derivative. 
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of equation (3.55) are geometrically conservative which means they satisfy the strong 
conservation form that we demand^'^. 



dt 



ipu'-ei 







By now we have established a mathematical framework to deal with nonlinear conser- 
vation laws in curvilinear coordinates. For practical astrophysical applications though 
we will have to dissect this framework a little further to time dependent coordinate sys- 
tems in order to work with adaptive grids as we motivated them for implicit numerical 
methods. 



3.2.3 Adaptive Grids 

In fluid dynamics we basically distinguish two main reference systems that suit unequally 
for various applications. The Eulerian frame is the fixed reference system of an external 
observer in which the fluid moves with velocity u whereas the Lagrangian specification 
describes the physics in the rest frame of the fluid. The transformation of an advection 
term for a density d between these two systems that obviously move with relative velocity 
u is given via the material derivative Dt = dtd + u • Vd. Hence, when we work with 
comoving frames, the coordinate system respectively the computational grid is time 
dependent. There is a number of purposes where strict Eulerian or Lagrangian grids are 
suboptimal and we need to generalize the concept of the comoving frame to ubiquitous 
grid movements. We will come back to this topic in chapter 4 when we deal with 
generation of moving grids in multiple dimensions. 

For now we want to evolve the strong conservation form for time dependent general 
coordinate systems. The time derivative of a density d in the coordinate system Sj-^) 
relative to a (e.g. static) coordinate system ^(a) is given by 9td(^) = dtd(^^^+V (^a)^ ^^-^{13) 
and from the view point of system S^^)) time derivative yields 

dtd^a) = 9id(^) - ±V(a)d. (3.56) 

With X we introduced the grid velocity and the second term on the right side we call 
grid advection, see e.g. Warsi (1981) [48] and Thompson, Warsi, Mastin [43]. Now let 
us consider an inhomogeneous advective term of a conservation law K = d^ -|- div (ud) 
referring to a fixed coordinate system. 

K = df - X • Vd -h div (ud) (3.57) 

When we apply transformation rules (3.50), (3.51), (3.52) to the spatial derivatives, we 
^°Here we simply used t^'ejBi ■ e' = f''ejSi' = t^^'ej. 
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gain the following form. 

K = dt - i 



\/ig|eM 








y^jg|"u • e'd 




vlg| 





(3.58) 



Admittedly this is not yet geometrically conservative, since it is not of an integral struc- 
ture. Following the idea of the Reynolds transport theorem we consider the temporal 
derivative of the volume respectively the determinant of the time dependent metric ten- 
sor Y^|g(x,t)| in order to study the conservation of a density function in variable volumes 
and obtain the strong conservation form for time dependent coordinate systems^^ . 



K = 9t J|g|d +d, 



|e* • (u — x)d 



(3.59) 



If we define the contravariant velocity components relative to the moving grid by C/* 
• (u — x) the upper equations yields 



g|K = at J|g|d +di A/|g|f/M . 



(3.60) 



3.2.4 Artificial Viscosity 

We introduced artificial viscosity in form of a three dimensional viscous pressure tensor 
in section 3.1.4 according to Tscharnuter and Winkler (1979) [ I !] by equation (3.20) but 
also mentioned that we will have to adapt it slightly for curvilinear coordinates. 

This has to do with one previous comment to the native form of tensorial quantities as 
they are generated in differential geometry. The crucial term in (3.20) is the symmetrized 
velocity gradient that accounts for shear stresses. The symmetrization rule is defined 
for lower indices and yields in component formulation [Vujg^^ = ^{ViUj + VjUi) which 
comes into conflict with the demand of vanishing trace TrQ = when the divergence 
term contains the unit tensor e. 

Proposition. The correct form of the viscous pressure tensor in general coordinates is 
Q = -'72^LcPmax(-divu,0) Vu ^ - ^gdivu^ . (3.61) 

Since this is an important result, we want to deduce this explicitly step by step. The 
viscous pressure tensor must be symmetrical by definition, i.e. Qij = Qji and as already 
stated TrQ = Q^^ = so we first consider its native covariant components 

Qij = -/i2max(-divu,0) (^{ViUj + VjUi) - ^gijdivuj (3.62) 



derivation of the strong conservative form for adaptive grids can be found in Appendix 6.1.4 
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where we renamed gi^visc/^ = M2- In order to proof that the trace vanishes, we need to 
raise an index with the metric and use the essential identity g''^gij = g^j = S^j. 

Q'j = Qij9^' = -Ai2max(-divu,0) Q/(V«Uj + VjU;) - ^5*jdivu 

The Ricci Lemma (3.36) for the fundamental tensor naturally also holds for the con- 
travariant components and we can permute g into the derivatives Vig^'^uj and Vjg^'^ui 
which yields twice the Divergence VjU* = divu when we conduct the contraction j ^ i. 
In three dimensions the summation 5% = 3 and we obtain our desired result 

Q\ = ... Q (2 divu) - ^3divu^ = 0. (3.63) 

Tscharnuter and Winklers form of Q is not compatible with our differential geometric 
definitions, since the symmetrization is only defined for lower indices whereas the unit 
element e of a metric space is only defined for mixed indices, meaning there is no such 
thing as 6ij. However, they and other authors like Mihalas and Mihalas [32] did not notice 
that little inconsistency since they considered mixed indices from the start respectively 
ID flows at last. Corrections have been suggested by Benson (1993) [■'>] albeit they do not 
epxlicitly concern curvilinear coordinates and were derived based on a TVD approach. 
The usefulness of our correction might not show until even more general coordinate 
systems respectively grids are considered. In chapter 4 we will discuss nonorthogonal 
grids and concern ourselves with the question at what step in the calculation geometric 
operations like raising or lowering indices should be made ideally. 

In the upcoming section we will recapitulate our general results for spherical coordi- 
nates without explicitly performing the discretization and marking down numerical fluxes 
according to (3.11) respectively Dorfi et al. (2004) [12]. It will become comprehensible 
why we suggest generating these discretized equations with Mathematica as well. 

3.3 Strong Conservation Form of RHD in Spherical and Polar 
Coordinates 

When we think of problem oriented coordinate systems in astrophysical applications 
spherical and polar coordinates are a first natural choice. In the following section we 
present the geometrically conservative form of the equations of radiation hydrodynamics 
in these cases exemplarily. We consider the transformation {x,y,z) — >■ (r E R"'",'!? E 
[0,7r],(^ E [0,27r]) with x = r sin -i? cos 99, y = r sin?? sin and z = rcos-d as spherical 
coordinates. The map {x,z) — ^ (r E M^,"!? E [0,27r]) shall be our polar coordinates 
X = r sin and z = r cos "!?. This convention was chosen in favor of somewhat simplified 
generation of the equations components for which we used the computer algebra system 
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Mathematica. Some of the equation generating files as well as some computational tests 
on consistency like Q*j = and Qij = Qji can be found in the Appendix 6.2. 

For the sake of completeness we mention some important geometrical attributes of our 
choice. 

Spherical Coordinates The covariant components of the metric tensor for these spher- 
ical coordinates are gij = diag(l, r^, sin i?) and due to their orthogonality the 
contravariant components simply yield g"^^ = diag(l, l/(r^ sin??)) and the 
volume element = sin t?. For the contravariant base vectors we obtain 



e*^ = (cos (/9 sin -i?, sin sin t?, cos 1?) 



,1? 



cos if cos I? cos 'd sin ip sin 

^ !y ty 

CSC sin (/9 cosy? CSC!? \ 

, ,0 (3.64) 

r r I 



and the non vanishing Christoffel symbols 



r^'tf^i = -r, T^'ip^ = -r SIT? i!) 

V^^r^ = V'P^r = \ r'^^^ = cot (3.65) 

Polar Coordinates Also polar coordinates are orthogonal and the metric components 
yield gij = diag(l,r^) respectively g^^ = diag(l, 1/r^). For the volume element we 
obtain \/|g[ = ^- The contravariant base vectors yield 



(sin T?, cos •&) 

( 

and the non vanishing Christoffel symbols are 



e'' = (3.66) 



T-d — r'^^r — r^^^ — — r. (3.67) 
r r 



3.3.1 Continuity Equation 

The strong conservation form of the conservation of mass (2.2) in form of the continuity 
equation 

5tp + div(pu) = (3.68) 
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for time dependend coordinates yields 



dt 



\P 



+ di 



|e* • (u - ±)p 



0. 



(3.69) 



Since the conservation of mass is scalar the divergence term may also be written in 
components S^v^jglp] + <9j[-\/|g|(u* — x*)/o] = 0. The explicit form of (3.69) in spherical 
coordinates yields 



d. 



sin dp 



+ dr 



sin [u!' — x^) p 



+ 



sin -d (^u^ — x"^^ p + d^p sin -d [u^ — x^) p 



(3.70) 



which we can rewrite via the transformation cos -d = p with differential d'd = — \/l — p^dp 
and divide the equation by \/l — p^ in the case of solely radial adaptivity. Note that 
even for solely radially adaptive grids the grid velocity x contains angular dependencies, 
e.g. i'' = fsini?cos(^ which obviously means x ^ f. 



dt 



r'^p 



{u^ — x^) p 



1 - p^ (u'^ - x^) p] + r^d^ [ (n^ - x^) p 







(3.71) 



With the radial grid velocity expanded we obtain the following form of the conservation 
of mass. 



r'^p 



r'^ iu^ — rJ 1 — p^ cos (/? ) p 



1 — p^ [u^ — rxl 1 — p'^ sin ip ) p 



{u'P-rp)p = (3.72) 



In the case of 2D polar coordinates the continuity equation yields 
dt 



rp 



+ dr 



r {vJ' — if) p — ry 1 — p^df^ 
respectively in case of radially time dependent coordinates 

dt rp + dr 



(u''-x']p 



r [u — r\ 1 — p^ ] p 



ryl - p^d^j, {u^ - fp^ p 



0. 



(3.73) 



(3.74) 



3.3.2 Equation of Motion 

The equation of motion that we consider in radiation hydrodynamics contains the con- 
servation of moment of plain fluid dynamics (2.4) as well as the radiative flux as coupling 
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term (2.18), gravitational force as indicated in section 2.3.2 and artificial viscosity (3.61). 



dt{pu) + div {puu + P + Q) + pG HRpH = 



(3.75) 



We investigate the elements of the energy stress tensor a little closer before we desig- 
nate the consistent strong conservation form. We define an effective tensor of gaseous 
momentum R that accounts for the motion of the coordinates 



R = r^^BiBj = p(u — x)u 



(3.76) 



with components r^^ = p{u^ — x^)uK The isotropic gas pressure tensor P as anticipated 
is defined by the scalar pressure and the metric tensor 



(3.77) 



with obvious decomposition P^^ = pg^^ . The viscous pressure tensor Q needs some 
special attention. We emanate from the following definition 



Q 



Ql'visc pes + 92^viscPmax(-div u, 0) 









( 


Vu 











(3.78) 



and remember that in order to get its contravariant components we need to raise them 
with our metric. 



Qij 
li 



(/ii -1-/22 max(-divu,0)) ( {ViUj +VjUi) - gij- div u 



Q'j = Qlj9' 



...[ g''{ViUj + VjUi) - S'j-divu 



■ ■■{ g'^g'^H^iUm + V,n„) - 5*^^divu 



1 



(3.79) 



We also promised to particularly calculate the trace of Q in spherical coordinates in 
order to confirm consistency. To do so, we need the mixed components which are listed 
on the following page. 
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Q r — Qrr9 ~^ Q'&r9 "I" QiprQ'^ — Qrr 

= — -div U + 9,Ur 
o 

Q^'& = Qr^g" + Qm9^^ + Qif'&g'''^ = Q'&r 

r 2 2 

Q ip = Qrifd + Q-d(pg + Qipip9^ — Qrip 
_ _l ^1. 

r 2 ^2 

V r 2 2 , 

= — ^div u+-Ur + \d^u^ 
6 J. 



Q^r ~ Qrr 9 ^ Q-drd ^ ~l~ QwrQ^^ ~ ^9 ^ o^Q(pr 



sm 



m^v \ r ^ 2 ^ 2 J 



sm 



cot + )^d^u^ + ^a^n^) 



1 

— ^div U + -Ur- + ^ cot 1?^^ + n ^'p'^'P 
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TrQ = (5\ = — -div u + drUr — -div u + -Ur + \d^u^ — -div u + -Ur + 
3 3 r 3 r 

1 1 

H K cot -i^li^ + o ■ n „ dMr 



1 



2 1 11 
—div u H 1 — ^ cot i9u^ + drUr H — ^d^u^ + 



-div u + ^ + cot i9u^ + a^n"^ + d^u^ + a^u'^ = 
r 

V ' 

div u 



9 . 2 Q ^'P^'P 



So in fact the viscous pressure tensors trace vanishes in spherical coordinates explictly^^. 

The native decompositions of our tensors R and P are contravariant hence we will 
consider also Q's contravariant components for the sake of coherency. In the following 
array we list Q^^ for spherical coordiantes. 



— -divU + drUr 

\ r 1 1 



1 



1 ^1. 
r2sin2^ V 2^"""^ 2 """^ 

r 2 2 



1 



-divu H — Ur -\ — T^d^u^ 



^ cot -Bu^ + ]^d^Ui) + ]^di)U^ 



r4 


sm^ 


1? 




1 






. 2 

sm 


1? 




1 






. 2 

sm 


1? 




1 






• 2 

sm^ 


1? 



1 



1 



1 



^- cot -i?n<^ + ^9<^u^ + ^<9tfn<^^ 
-divu H — H — K cot i7ti^ + 



9 . 2 Q ^•fi^'p 

sm^ if ^ ^ 



We will investigate the following strong conservation form of the equation of motion 



d. 



\pu 



+ d, 

pdi v/|g|^>e 



g|(R + P + Q) -e 

4tt 



+ 

K/JA/lglpH 



0. 



(3.80) 



The fc-th component of the strong conservation equation of motion is given by the /c-th 
cartesic component of our unit vector, e.g. e*" = cos sin t?e^ -|-sin (p sin i?e^ -|-cos f?e^ and 



^^The geometrical terms 2/r and coti9 are Christoffel symbols of the spherical coordinates (3.65) from 
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its derivatives. The projection of each physical tensor on the contravariant coordinate 
hnes can be simphfied by its contravariant components with respect to its covariant basis 
without losing strong conservation form, i.e. f • e* = /^e^. We prefer this form with 
contravariant components since it meets the native design of the stress tensor (2.3). 



{a. 



\pu en 



^in ^ pin ^ Qin\ ^ 



pdn Vlgl^'e" 







(3.81) 



The three components of upper equation in spherical coordinates are listed in the fol- 
lowing arrays. 



dt pr^ 1 — fi^ cos ipu^ + pr^ /ip cos ipu^ — pr^ ^ \ — p?- sin ipu^ 

-r-='yi^^sinv3(^p«'^(M''-i'') + Q'''^^ - 9^ r'^{\ ~ p^) cos ip(ypu'' {u^ - x^) + Q^''^^ 
:^p^\~p^ cos Lp[pu'>{u^ - i'') + -ip + Q''''^ - r\\ - p'')sui^[pu^(u^ - x^) + Q^^^ + 
'^J\- p^ cos ip[pu''{u'P - i*") + 0'^''^ ^r'^p cos ^[pu^{uf - i*' + Q'^'')^- 
-r^ -v/r^-^sin ( PU^(m'^ - x^) + ^ 5tP + Q'' 



+p^9r 1''^ \J ~ P^ COS ip^ ~ pr^l — p^ cos tp^ 



r sin (/3<I> 



An 



^y/ 1 — p^ cos ipH^ + pr cos ipH^ — \J \ — p^r sin cpH'^^ 



0. 

(3.82) 



dt 



-\/ 1 — p'^p sin ipu^ + r^pp sin ipu^ + r'^ -\/ 1 — p'^p cos (fu^ 
(^p«'-K - i'-) + P + +r''psin^{^pu^{u' - i'-) + Q^^^ + 

(1 - At2) sin v5(^p«'-{«'' - i'') + Q-^"^ + 



^ \J \ — p^ sin (/3 
+r3v^l - p2 cosi^(^pu^{m'' - i'') + Q'''^^ 
+rVv/l - sin¥:.(^p«''{«'' - + ^ + Q''''^ + r3(l - ^2) cos (^pu*" - i'') + Q''^^ 
^A/r^sin</p(^p«'X«'^ - x'^) + Q"^''^ + rVsin</3(^p«''{M'^ - x"^) + Q'^''^ 
-^/l — /i^ cos (p ^pii*^ (m''' — i''') + 



1 — p2 sin (^(j> 



/ir — ■ p2 sin (p<I> 



r-2(l - p2) 
1 



p- 



r cos ip<E> 



4-71 



r2 ^ y/ 1 — p^ sin ip//' + pr sin (p_ff + \/ 1 — p^r cos (fH'^^ 



0. 

(3.83) 
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+9, 



r(l - /i^)* - (^f.H'- - r^/l^H^^ 



+dr 



(3.84) 



When we paste the contravariant components of the artificial viscosity exphcitly, the 
1-component yields as an example 



dt 



pr^(l — /i^) cos (fiu^ + pr''' ^yj 1 — fi^ p cos ^pvf — pr^ (1 — p^) sin ipu^ 



+dr 



^(1 — p,^) COS ip(^pu^ {u^ — x^) + p — -div U + drUr 
p,\/ 1 — p? COS 1 pu'' {u^ — X^) H ^ ( tltf H drU^ H (?^Mr 

\ \ r 2 2 



r^(l — p-^) \ r 2 2 



— -v/l — /i^(9p r^(l — /^^) cos ip( pu'"(m'' — X-'') H — ^ ( H — drU^ ^ — di)Ur 

L \ V r 2 2 

+r^p\/l — p^ cos ifil pu^{u^ - X-'') + \p + \ ( — -divu + -Ur + \difU.ff 
\ V 3 r r-^ 

-r^{l - n'^)sin(p(puf{u''' - d;'') + — — !^ — 5- f , ^ My + ^9^'"'? + 7;9a'^v>) 

\ r4(l-/i^)V y'i_^2 2 2 y 

r^(l - p'^) cos ip ( pu^{u'^ — i*^) + — ( — -Uu, + -dwUr + -drU^ 

L V r^(l — /i^j V r 2 2 

-r^/i-y/l - p'^cosip(pu^{u'f - xV) + — — — — ( , ^ My + + 



— (1 — p^) sin j pu'^ (u'^ — i''' 



r2(l-/i2) 



P+ 



1/1,. 1 1 M 1 ^ 
— ^ - -divu+ -u^ + — + — ^dv-^v 



^(1 - /i2)coS(/3$ 



A^/T^-Ti^ 9(j /^r ^ \ — p^ cos 



- r sin 



47r 



)- 

— jU^^-^l — /i2 cos ifH'^ + pr cos^pH^ — \J ^ — p^r sin ipH'^^ 



(3.85) 



The two components of the equation of motion in polar coordinates yield 
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+p(^y/l^^dr r* - - M^Sm ) - ^r(^^Jl- p?H'- + rpH 



(3.86) 



a* 



r 2 

r/^p« — r 



2 „,,!» 



+p(^p9, r* + ^1 - p2a^ ^1 - p2$ ^ - ^r{pH'- - r^l - ii'^H^^ 



(3.87) 



3.3.3 Equation of Internal Energy 



The energy equation we are going to impose contains the thermodynamics of the fluid 
(2.5) as well as the energy exchange with the radiation field (2.18) and viscous energy 
dissipation, expressed by the contraction of the viscosity with the velocity gradient Q : 
Vu. 

dt{pe) + div (upe) + P : Vu - 4iTKpp{J - S") + Q : Vu = (3.88) 

Since we assume isotropic gas pressure P = pg we can reformulate its contribution via 
the Ricci Lemma (3.36) and obtain a very simple scalar expression. 

P : Vu = g'^pViUj 

= pViu' =pdiYU (3.89) 

The viscous energy dissipation is a little more work to do, since both matrices contain 
a priori no zeros. However, for this contraction we do not need to appeal to strong 
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conservation differentiation since tliis term does not fiave its seeds in a conservation law. 



jUk 



Q : Vu = Q'^V.Uj 

= Q"'drUr + Q''^(drU.i) - ^U-d^ + Q""^ [drU^ - + ... 

= Ediss (3.90) 



Of course also this quantity from now on denominated E^iss is calculated with the aid of 
Mathematica. 



The strong conservative form of the energy equation yields 



|g|/5ee* • (u — x) +Y|g|pdivu — 
■At:J\%\kpp{J - S) + A/lglEdiss 



0. 



(3.91) 



We project the velocities on the covariant base, multiply with \/jg| and divide by \/l — /x^ 
to obtain the energy equation in spherical 



dt 



r^pe 



1 - p^peiu"^ - x"^) 



+ 



pe{u^ — x^) + r^pdivu — 47rr^Kpp(J — S) + r^Ejiss 



and polar coordinates. 



dt 



rpe 



rpe{vJ~ — xf) — rJ 1 — p?dfj_ pe{u — x ) 



+ 



(3.92) 



rpdiv u — 47rrKpp( J — S) + rEdiss = (3.93) 



3.3.4 Equation of Radiation Energy 



The simplified frequency integrated radiation energy equation in the comoving frame 
was defined in equation (2.28). 

atJ + div(uJ) +cdivH + K : Vu + cxp(J-5) = (3.94) 
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For the scalar energy input of radiative pressure into the material we analogously define 
a new coupling variable 



K : Vu 



coup 



In strong conservation form it yields 



dt 



\J 



e'' • ( J(u - x) + cH) 



+ 



|g|Pcoup + Vlg|cXp(-f-5) =0. 

In spherical coordinates with fixed "Q we obtain 



dt 



+ dr 



r du 



+ 



r 



and the radiation energy equation in polar coordinates yields 



dt 



rJ 



r(JK-x^) + cF'^) 
+ r-Pcoup + rcxp{J - S) = 0. 



(3.95) 



(3.96) 



(3.97) 



(3.98) 



3.3.5 Radiative Flux Equation 

The radiative flux equation, another vectorial conservation law remains to be sketched in 
strong conservation form for general non-steady coordinates. We emanate from equation 
(2.28). 

dtU + div (uH) + c div K + H • Vu + cxrH = (3.99) 

and for a start define an effective radiative flux tensor L analogously to the effective 
tensor of gaseous momentum. 



(u -x)H 

= {u'-x')W 



(3.100) 
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The contribution of radiative momentum to the material H • Vu will be abbreviated via 
another new coupling variable with components 



The geometrically conservative form in non-steady coordinates of upper equation will 
be considered in the following shape. 



The three components of the radiative flux equation in spherical coordinates yields 




(3.101) 
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dt 



„3 



1/1 -/i2 COS <p(^H'-{u'^ -*'■) + c/f''^ 



COS i^if + cos v'i?'' - - //2 sin (^if"^ 



- r^/icos ipl 



+ r^lJ.\/ 1 — fi^ cos </3 



cos - i"^) + cK'f''^ + cos {h^ {uf - x^) + cK'^^^ 

sin ^ (hv (w"^ -x'^)+ cKf"^^ 



cos v'F^oup + r/^ cos ipV^o^^ - ryJT- fi'^ sin tpFtonp ] + 



+ constr^ I Y 1 



0. 



(3.102) 



Vl -M^ sin(^^ff(u'' - i'') + ci^''''^ + rVsin<^(^/f''(u'' - i'' ) + cE"''' ^ + 
f r V 1 - cos (^i? (x'' - ) + cK''"^ j -9^ r2(l - m^) sin¥'^i?'-(u'' - i'') + 0/^'''"^ + 
+ r^^l^/l- H'^ sin ip(h^{u^ - x^) + cK^-^^ +r^{l- n"^) cos tp^m^u^ - x'^) + cK'^f^ + 
^ -/T-^ sin <p ^H'-iuf -xV) + cK'f^^ + r^n sin ip (h^ (uf - x'^) + cK'f^^ + 



^ {^\J I - IJ.'^ sin ip^laxtp + A"" sin V'Fcoup +r\/l- iJ?cos ^^ton^ + 
+constr^ ^-s/ 1 — /i^ sin ipH'^ + ^ \ — \x'^ cos i^if''' + rjx sin (pH^^ 



0. 

(3.103) 



0. (3.104) 
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The two components in polar coordinates can be found in the following arrays. 



+ 



constr[ W 1 — iJ?H^ + r^iH^ 



0. (3.105) 



dt 



+ 



0. (3.106) 
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In the forthcoming sections we want to motivate and discuss the necessity for multi- 
dimensionally adaptive grids for astrophysical apphcations with radiation hydrodynamics 
in 2D and 3D and thereby justify our excessive mathematics hitherto. At first we will 
give a brief introduction to fundamental concepts in grid generation and then phrase 
three crucial postulations for suitable adaptive grids in implicit RHD numerics. We will 
discover that not all three postulation can be fulfilled and therefore occupy ourselves 
with possible compromises and prospective continuation of our ideas and approaches. 

4.1 Basic Principles of Grid Generation 

Finite difference and finite element methods for numerical solutions of partial differential 
equations require a grid which represents the physical domain in discrete terms, often 
referred to as the logical space. Generally the grid is a set of points, lines, (hyper)surfaces 
and (hyper)volumes that arise from an appropriate transformation between these two 
domains. In many applications the points, also called nodes and volumes or cells are 
to be chosen in a way that the underlying problem gets simplified; for instance a trans- 
formation to a rectangular, an orthogonal or a merely smooth domain. Such structured 
meshes can be obtained via coordinate transformations, algebraic methods like poly- 
nomial interpolation or differential methods respectively hybrid variants (a textbook 
giving a profound overview is e.g. Knupp and Steinberg [22]). In chapter 3 we antic- 
ipated generating a structured mesh by discretization of coordinates we obtained from 
a well known coordinate transformation even though we never explicitly conducted this 
transition from (r, 93) to (rj, i?^, and differential to difference operators. 

Unstructured grids cover the computational domain with cells of arbitrary shape and 
do rather appear with complex geometrical domains. Practical problems of unstructured 
meshes arise from numbering and ordering the cells and edges which puts high compu- 
tational demands even if boundaries do not move. Clearly also the numerical solution 
of a partial differential equation on such an unstructured is more expensive than on a 
structured mesh. However, one advantage of unstructured grids is a quite handy local 
mesh refinement by dividing grid cells appropriately (Liseikin [30]). 
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Grid generation itself must be understood as an essential step of finding suitable 
methods for solving PDEs numerically^. In the following sections we want to point out 
some possible approaches to grid generation in 2 and 3 dimensions. However, for the 
sake of clearness 2D grids will be favored. 

Time dependent problems with steep gradients as they appear in hydrodynamics state 
special challenges for numerical techniques, likewise for numerical grid generation. Ba- 
sically one can distinguish between three main approaches to handle non static finite 
differences or finite elements; the Eulerian, the Lagrangian and methods that allow 
generally adaptive grids. As pointed out in section 3.2.3, Eulerian grids remain static 
throughout the computation and all physical quantities and fluxes are calculated with 
respect to steady nodes and cells whereas Lagrangian grids move with the very velocity of 
one particular physically interesting component, in most cases the gas. Analysis reveales 
(Furzeland et al. (1990) [15]) that generally adaptive grids are definitely privileged. 
Connecting the grid with physical quantities respectively some measure of their compu- 
tational error and adjusting its profile iteratively provides the best results whenever the 
number of nodes is fixed (i.e. no mesh refinement). 

Fundamental ideas of dynamically adaptive grids in RHD were developed by Winkler 
(1977) ], Winkler, Norman and Mihalas (1984) [ ] and [33] respectively a sophis- 
ticated method for grid control in ID was found by Dorfi and Drury (1987) [II]. In 
the following paragraphs we shall investigate some potential amplifications to multi- 
dimensionally adaptive grids. 

4.2 The Grid We Are Looking For 

First and foremost we are looking for a grid that geometrically operates radiation hydro- 
dynamics in problem-oriented domains. Since we know that all stars rotate and therefore 
are rather geoids than spheres, polar and spherical coordinates do not match optimally. 
In latter geometries, shock waves would not propagate alongside coordinate lines respec- 
tively within coordinate surfaces but skew to them. Immense fluxes in multiple directions 
would be the consequence which is numerically troublesome. Hence, an adaptive grid 
that satisfies our demands ideally would have to allow adaptiveness in more than one 
dimension but remain orthogonal all along. 

Postulation. The grid should be adaptive in more than one dimension but remain or- 
thogonal throughout the computation. 

In the case of strong conservation finite volume methods in complex geometries the 
orthogonality of the coordinate lines would ensure that the metric tensor (as we have 

^ There are popular numerical methods working without meshes like the family of spectral methods or 
smoothed-particle hydrodnamics but as the chapter title suggests we concern ourselves with grid- 
based techniques exclusively. 



64 



4.3 Remarks on Grid Generation Methods 



seen a function of space and time) remains diagonal. With more general metrics, not only 
the volume element |g| becomes rather elaborate as it is calculated by the determinant 
of the metric. All projections and summations of our physical tensors require far more 
arithmetic operations. For instance, each of the nine contravariant component of the 
viscous pressure Q requires nine arithmetic operations. In discrete terms this implicates 
to execute these 27 x 2 operations at each node and time step at worst. 



With orthogonal grids, all off diagonal components of the metric are zero g^^ = i ^ j 
and in this case 54 operations reduce to effective six. After all this objection is not 
that relevant practically, since in implicit radiation hydrodynamics most CPU time is 
consumed by the matrix inversion (Dorfi (1999) [ ]). 

Even if orthogonality can not be maintained (which indeed will be one of our con- 
clusions), the adaptive grids by any means needs some kind of mollifying property that 
avoids excessively twisted, acute-angled or one another overtaking cells. 

Postulation. Grid cells should basically maintain their local geometries. 

Moreover, the one dimensional special case of simple radial geometry shall be inherent 
in more dimensional grid generation so that one can easily compare calculations, e.g. 
with previous works in implicit RHD numerics like Dorfi (1999) [I '-!] respectively Dorfi 
et al. (2004) [ ]. 

Postulation. The more dimensional grid should contain the one dimensional special 
case of radial geometry. 

In the upcoming section we will discuss several approaches to grid generation in the face 
of problem oriented geometries and suggest feasible approaches to multi-dimensionally 
adaptive grids in 2D and 3D radiation hydrodynamics. 

4.3 Remarks on Grid Generation Methods 

In this section we present some possible approaches to grid generation for implicit adap- 
tive mesh techniques and discuss their compatibility with our three postulations. A 
common premise is an a priori fixed number of grid points that is redistributed during 
the temporal evolution. Again our perspective emphasizes on methods for neither fully 
Eulerian nor fully Lagrangian grid but a hybrid where gas velocity and grid velocity are 
associated but not equal as introduced in section 3.2.3. 

Before we glance at adaptive grid techniques we present some fundamental static grid 
generation methods that could be adapted for time dependent problems. 




H mj 
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4 Orthogonal and Nonorthogonal Adaptive Grids 

4.3.1 DifFerential Methods 

Grid generation using differential equation tecliniques are very popular for structured 
meshes, especially with complex geometries. The idea is to set up a system of partial 
differential equations whose solution defines a coordinate transformation. Well known 
properties of special differential equations can be used to control coordinates to some 
extent and on the other hand help to find suitable solution methods. These and a lot of 
other methods are well described in several text books, e.g. Castillo [8] or Liseikin [30]. 

Grid generation methods based on solving partial differential equations are natu- 
rally classified elliptic, hyperbolic and parabolic. Hyperbolic equations are preferred 
for meshes around bodies but prove inappropriate for internal and closed domains which 
means they naturally disqualify for astrophysical applications. Moreover, inherent sin- 
gularities propagate and grid oscillation and overlapping of volumes is encountered fre- 
quently. Parabolic techniques combine some positive aspects of hyperbolic properties 
like being able to formulate an initial value problem on the one hand and elliptical prop- 
erties like diffusion effects smoothing out singularities on the other hand. Numerous 
references to hyperbolic and parabolic as well as combined techniques can also be found 
in [30], pl91 ff. Mathematical fundamentals as well as research activities in the field of 
differential grid generation are presented in Knupp and Steinberg [■^■^]. 

Exemplary we want to sketch some basic properties of elliptic systems using Laplace or 
Poisson equations. Elliptic equations have one important characteristic as they ensure 
smooth solutions in the total interior of the domain even if the boundaries are non- 
smooth. For smooth boundaries continuous derivatives of the resulting coordinates are 
also granted on the boundaries. Elliptic equations that obey the extremum principle like 
the Laplace Equation are also known to have very small tendency for folding of grid cells 
(see e.g. Spekreijse (1995) [ ; ], Zhang et al. (2006) [53]). However, we should mention 
that elliptic systems are numerically costly and Laplace-type equations allow practically 
no control of the geometric properties of the grid. In anticipation of more sophisticated 
methods the Laplace system is presented briefly. 

4.3.2 Laplace Systems 

Many of the commonly used techniques to generate grids via differential equations are 
derivations of a method proposed by Winslow (1966) [52] where the nonlinear two di- 
mensional Poisson equation is solved. In these approaches the constructed meshes can 
be understood as equipotential surfaces. 

We shall study an archetype of such an elliptical grid generation method and consider 
a computational domain H*^, a physical domain X'^ and a coordinate system 

x(0^(x(,/(6), i = l,.--,d (4.1) 
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4.3 Remarks on Grid Generation Methods 

The most obvious elliptic system is the Laplace system that can be formulated in the 
physical 

Ax' = r.x' = (4.2) 

or computational domain (interchange x and ^) which is equivalent with g''"9/5mx' = 0, 
a coupled quasilinear system of partial differential equations in H. In most cases the 
transformation x(^) will be considered since this formulation in the physical space allows 
direct control of grid spacing and orthogonality. 

One main requirement for the transformation is that the Jacobian of the transfor- 
mation x(^) : — >■ is nonzero, i.e. the two domains are diffeomorphic which is 
confirmed by the theorem of Rado for harmonic functions. This leads us to the first 
discrepancies with our postulations made before. There is no diffeomorphism between 
Cartesian and polar respectively spherical coordinates that includes r = 0. Moreover, 
even if we cut out the origin, polar and spherical coordinates do not even satisfy the 
Laplace nor a Poisson equation. 

Nevertheless we briefly discuss one application and consider the transformation x(^, rf) 
in a a two dimensional simply connected bounded domain and let r? be a unit square 
as commonly imposed. In section 3.2.1 when we introduced basic differential geometric 
relations, we also mentioned that the metric tensor g is symmetric which leads to the 
following boundary value problem^. 

(/«a| + 2g^^di,dr, + g^^d'^^x' = (4.3) 

Computations as well as differential geometric deliberations reveal that grid lines ob- 
tained by (4.3) are attracted to concave boundaries and repelled in the vicinity of con- 
vex boundaries. In stellar astrophysics this would implicate that the mesh gets tighter 
towards the center of the star automatically. The bottom line is that such elliptical 
systems as pure boundary value problems disqualify for grid adaptiveness with respect 
to inner domain dynamics. 

4.3.3 Conformal Mapping 

When thinking about our first postulation of more dimensional adaptivity and orthog- 
onality one might consider conformal maps since they are angle-preserving. Conformal 
transformations would remain orthogonal throughout the evolution of the grid which is 
represented by a diagonal metric gij = diag{gii). Unfortunately coordinate caustics and 
overlapping problems are not the only challenges one would have to face. Both domains, 

^For our orthogonal reference, the polar coordinates we obtain a contradiction (9^ + 
ia|)(rcosi9,rsini9) / 0. 
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4 Orthogonal and Nonorthogonal Adaptive Grids 

and X" have to be conformally equivalent which is expressed by the fact that the 
metric tensor is a multiple of the unit matrix. 

9v = 9^3 (4-4) 

Again we come into conflict with our postulations given that our reference systems, polar 
and spherical coordinates are not conformally equivalent to Cartesian ones^. 

4.3.4 Algebraic Methods 

Algebraic grid generation methods rely on coordinate transformations constructed by 
interpolation techniques. According to concrete requirements to the solutions, different 
interpolation functions respectively polynomials are evaluated and transformations are 
governed by the coefficients of these polynomials. 

These schemes prove computationally efficient but lack some convenient intrinsic prop- 
erties we presented in context of differential methods. Algebraic grid generation requires 
significant control techniques to obtain practicable meshes otherwise skewness of cells 
and grid point spacing can get messy. There are several ideas to regulate the grid such 
as adding differential calculus elements in order to inhibit that coordinate surfaces come 
too close or even overlap respectively degenerate. Other approaches directly influence 
the interpolation coefficients. 

In the following section we present an idea of a hybrid method that combines concepts 
of algebraic methods with equidistribution techniques. 

4.3.5 Equidistribution 

The idea of equidistribution methods is to place the grid nodes in a way that the numerical 
error is distributed uniformly throughout the computational domain. In practice that 
means to adjust the mesh size locally according to the local numerical derivatives - the 
steeper the gradients the higher the spatial resolution required in order to get uniform 
resolution. In one dimension a given physical quantity gets equidistributed along an arc 
(see e.g. Thompson, Warsi, Mastin [1'^], Dorfi and Drury (1987) [ ]) and the grid is 
determined uniquely by imposing a proper grid equation. 

The approach presented in [1 i] is to construct a grid equation distributing the nodes ui 
proportionally to a desired resolution R which is defined in terms of the arc- length along 
the graph of an approximated function f, R = \/l + (df /dx)"^ which originates from the 
definition of the line element ds^ = gijdx^dx^ . The choice of function f depends on what 
regime in the calculation shall be the focus of resolution, e.g. pressure or temperature. 

^In polar coordinates the metric yields Qij = diag(l,r) / g{r,d)5^j. 
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4.4 Equidistribution and Interpolation Attempt 

A smoothing respectively variation limiting control of the grid is obtained by imposing 
constraints on grid spacing and the resolution function. 

Generalizations to more dimensions follow different strategies. Several attempts have 
been made to generalize one-dimensional equidistribution directly to two and more di- 
mensions by applying a series of ID equidistributions along coordinate arcs. The concept 
of arc equidistribution however can only be satisfied locally but not globally in the do- 
main as pointed out by Anderson (1987) [ ] as well as Hunag and Sloan (1995) [20]. 
Moreover, successive arc equidistributions along fixed coordinates fails to adapt the ge- 
ometries of the cells which was one of our postulations. Hence such a straight approach 
is rather comparable to adaptive mesh refinement methods. 

As motivated in previous chapters, it is reasonable to perceive physical quantities and 
their numerical pendants in strong conservation form in a volume integrated framework 
respectively even in terms of differential forms. So it would seem natural to investigate 
methods that equidistribute a quantity with respect to surfaces and volumes, studied 
e.g. in Delzanno et al. (2008) ["] or Sulman, Williams and Russell (2009) [41]. 

One of our first attempts to obtain a multi-dimensionally adaptive grid that remains 
orthogonal during temporal evolution was to combine one dimensional equidistribution 
techniques and algebraic methods. 

4.4 Equidistribution and Interpolation Attempt 

The intention is to generate a 2D grid that allows certain deformations from a polar 
coordinate system but remain orthogonal globally^. One family of curves shall represent 
the angular coordinate and a family of normal curves would define the radial one. Op- 
timally we would start with a polar grid and according to an equidistribution principle, 
e.g. the angular curves would deform and intrinsically define the family of orthogonal 
radial coordinates. 

4.4.1 Toy Model 

Let us consider a transformation (^(x, y, t), r]{x, y, t) where the coordinates (^, rf) are rep- 
resented as angular and a radial family of curves. In favor of simplification we assume 
azimuthal symmetry and confine our analysis on a quadrant of the real plane. The most 
obvious deformation we would want to be considered is oblateness by rotation which 
means in 2D the deformation from a circle to an elliptical shape. As naive toy model we 

*For the time being we shall blind out that this undertaking is futile altogether and refer to posterior 
findings. 
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4 Orthogonal and Nonorthogonal Adaptive Grids 




Figure 4.1: Angular coordinate lines of toy model 




Figure 4.2: Coverage and Caustic 

let ?y : M X M X — M be a polynomial with a family of time dependent coefficients^. 

v{x,T,t) = J2 ai{T,t)x' (4.5) 

The coefficients a, are partially determined by boundary conditions as we demand normal 
tangents of r/ at y = and a; = in order to guarantee smooth transitions to the other 
quadrants of the real plane. The ^ curve at a point xq is expectedly determined by 

C(x,r,t) = r/(xo,r,t) — -{x-xq). (4.6) 

r/'(xo,r, t) 

For principal analysis we set imax = 3 and studied this system in Mathematica. This 
choice determines the angular curves with noted boundary conditions uniquely, each 
further term would bring ambiguousness into the grid that could be used to include 
physics respectively grid control mechanisms. 



^Evidently with this ansatz our reference frame, the polar coordinates are not included which we also 
blind out for the present. 
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4.4 Equidistribution and Interpolation Attempt 




Figure 4.3: Normal vectors 



4.4.2 Analysis 

We define differently equidistributed points at the x- and the y-axis. Since the boundary 
condition for the rj curves at y = would imply infinite slope, we rotated the whole 
setting by a = 7r/4. What we get is a set of curves that increases in ellipticity with the 
radius which describes the physical situation of a rotating configuration, drafted in figure 
4.1. As already mentioned both the angular coordinate lines and their normal curves 
are determined a priori. With determination of one orthogonal curve at a point xq, it 
necessarily would have to be normal at all intersections if there existed such a family 
of curves. The most obvious observation we make in figure 4.2 is that the radial curves 
do not intersect with the natural point of origin nor reproduce the (45 degree rotated) 
X-axis and cause coordinate caustics. Of course with growing order of the polynomials 
further boundary conditions can be imposed. 

However, this leads to the second quite evident observation that the orthogonal curves 
to the outermost polynomial are not normal at all other intersections with inner polyno- 
mials. If we wanted to impose orthogonality at n intersections, the order of polynomial 
in this ansatz would at least have to match the number of radial nodes imax ^ n which 
is rather pointless when we consider e.g. 500 radii. 

This graphical analysis already discloses the limitation of grid generation in two di- 
mensions considering our postulations. Any deviation from polar coordinates reproduces 
either nonorthogonal grids or fails to reproduce the domains boundaries; hence postu- 
lation one will have to be renounced in its original form. The only kind of orthogonal 
adaptiveness in a polar grid is dislocation of nodes alongside the original coordinate lines 
which corresponds to adaptive mesh refinement methods. In the following section we 
want to substantiate our conclusions analytically. 
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4 Orthogonal and Nonorthogonal Adaptive Grids 

4.5 On the Peculiarity of Orthogonal Coordinates 



In section 3.2.1 we studied differential geometry of coordinate transformations such as 
the definition of base vectors (3.27) and the metric tensor (3.30). We noticed that orthog- 
onal coordinates imply diagonal metrics, a factor that can be used to study this family 
of coordinate transformations. Since one feature of coordinate lines is that all other 
coordinates are constant alongside the coordinate, transformations have to factorize and 
the following product ansatz is justified*'. 

x(e,r/) = a(Oa(r/), r/) = 6(e)/3(r/) (4.7) 

With aforesaid definitions we obtain the base vectors = (Qa',/36'), = {aa'bf3') and 
following covariant metric components. 

g^^ = a a + p b 

2 /2 I 1,2 a>2 

griT) = a Oi ^ + h ji 
9(ri — 9vi ~ flQ!a'a' + b(3b' [3' (4.8) 

In order to obtain a PDE we can study and solve analytically, the off diagonal entries 
g^jj are set zero and functions of are separated from functions of r] which yields 

^ = -^ = A (49) 
bb' aa' ^ ' 

where A G M. With the substitution A = and = 6^ we obtain A' / B' = A repectively 
0? = Xb^ + c with c G M is an arbitrary constant. Our boundary conditions come from 
the postulation that the coordinate system shall in principle contain the special case 
of polar coordinates hence we set the point of origin for the radial coordinate to zero 
a(0) = 0, 6(0) = and obtain 

a = \/A6 (4.10) 

which means that the a and b differ just by an arbitrary factor. When we replace 
= A, fi"^ = B we obtain B = —\A + k and determine the constant k with the 
boundary conditions that the angular coordinate reproduces the x-axis at r/ = and the 
y-axis at ry = 1 to A; = A yielding 

/3 = VaVi -a2. (4.11) 



'We confine the analysis on 2D but reasoning and results are equivalent in 3D of course. 
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4.6 Quasi-Polar Coordinates Attempt 

A general orthogonal coordinate transformation which contains the reference frame of 
polar coordinates hence necessarily has to satisfy 

r,) = ^HOaiv), r?) = ^b{0 ^ I - a{r,r . (4.12) 

The scaling factor ^/\ can evidently be set to ^/X = 1 without loss of generality which 
leads to our result 

a = 6, a^ + ^^^l. (4.13) 

The postulation for orthogonality inevitably leads to radial symmetry and the grid can 
never gain an elliptic or even more deformed shape. 



4.6 Quasi-Polar Coordinates Attempt 

Evidently it is necessary to adapt our postulations for a multi-dimensionally adaptive 
grid in the light of our previous results. Since orthogonality can not be sustained we 
renounce that strong demand and concentrate on producing adaptive grids that include 
our reference frame. The following ansatz includes algebraic elements as well as some 
differential boundary conditions that control the grid properly. 

4.6.1 Toy Model 

Again we consider a transformation ^{x,y,t),ri{x,y,t) with radial and angular coordi- 
nates (^, r)) and azimuthal symmetry. The product ansatz we suggest is 

^'max ?max 

a(e,t) = X^a,(t)e, b{^,t)=J2bi{t)e 

i=l 1=1 

Jmax J max 

a{ri,t)= (l+J2»jitW)smri, Piri,t) = {l + (tW) cos rj 

(4.14) 

where (a^, hi, aj, Pj) are time dependent coefficients that govern the adaptivity of the grid. 
One main feature of this approach is that it allows to turn off all deviations from radial 
symmetry by setting (ai,6i) = 1 and (ai,bi) = OVi > 2 plus (aj,l3j) = which yields 
ploar coordinates. Basically the unknown coefficients are determined via appropriate 
boundary conditions and whatever coupling with physics would be implemented (e.g. 
equidistribution along coordinate arcs or with respect to cell volumes). 

In order to discuss basic characterictics of this proposal, we set imax = 2 and jmax = 3 
and dissect the system in Mathematica. 
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(a) 



(b) 





(c) (d) 
Figure 4.4: Toy Model Grids 



4.6.2 Analysis 

For feasible analysis of course we need to constrain the degrees of freedom of our toy 
model but leave some undetermined parameters to be manipulated. With boundary 
conditions that on the one hand assure smooth transitions to other quadrants (normal 
derivatives at x- and y-axis) and on the other hand are consistent with the reference 
frame we effectively solve a system of algebraic equations in eight unknowns and gain 
three free parameters. The Mathematica Notebook we used for solving and plotting these 
functions can be found in Appendix 6.2. 

One example of a coordinate transformation that in principle satisfies all postulations 
from 4.2 besides orthogonality yields 

X = ^ cos T] 



^^'^ + ) V + 4vr(l + h,e) ^ + 



+ 



4vr(l + 62O 
-/33vr3 + 462C - b2/33vr3e) 2 



^2(1 + 626 



(4.15) 



where 61, 62 and /Ss G M are free parameters. Certainly they are not totally free but need 
to be confined to certain bounds that ensure unique covering of the domain. Further 
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4.7 Strong Conservation Form in Nonorthogonal Coordinates 



Table 4.1: Maximaa of geometrical parameters for toy model 





61 


62 


/33 


max|c/i2(C,r/)/^| 


max^^^j(^,?7) 


(a) 


1 











1 


(b) 


0.8 








0.18 at (0.85,0.73) 


0.64 at (0.9,0.71) 


(c) 


0.8 





0.5 


0.13 at (0.89,0.79) 


1.19 at (0.86,0.33) 


(d) 


0.8 


0.2 


0.5 


0.25 at (1,0.43) 


1.62 at (1,0.29) 



constraints emerge from differential geometric requirements to the grid like a maximal 
skewness of cells. 

In figure 4.4 we visualize grids with varying parameters that are listed in table 4.1 along 
with some characteristics of the grid. The off-diagonal metric elements gi2 = §21 7^ 
quantify the skewness of the grid and v^j,gj designates the local cell volume in rela- 
tion to polar coordinates. Both characteristics were numerically maximized in Math- 
ematica in the ranges G (0,1] and rj S [0,7r/2] since the functional expressions get 
rather bulky. For the simplest quasi-polar model (b) however the coordinate trans- 
formation is X = i^cosr], y = 0.8.^ sin 77 and the grid characteristics yield explicitly 
512, (b) = — 0.36^ cos 7? sin r/ and \/gj,gi = 0.64. Due to our boundary conditions the 
minima are mostly to be found at the origin respectively at the axes, for details see 
section 6.2. 

4.7 Strong Conservation Form in Nonorthogonal Coordinates 

In section 3.3 the equations of RHD were explicitly formulated in polar and spherical 
coordinates and we noticed that even with orthogonal coordinates the expressions get 
rather cumbersome. From the view point of tensor calculus, the number of operations 
like raising and lowering of indices with non-diagonal metrics increases with the square 
of the spatial dimension considered. Strong conservation form radiation hydrodynamics 
in general curvilinear coordinates hence can be expected to arouse extensive equations. 

In fact, if we take the full ansatz (4.14) into our Mathematica files that generate the 
equations of radiation hydrodynamics in strong conservation form the computing time for 
just generating the equations gets sizable and the outputs become huge. The question 
arises at what point of such calculations the geometric relations and terms are to be 
determined respectively computed ideally. A profound answer will have to be given in 
future investigations. However, without a doubt the actual solution of the RHD system 
which is a matrix inversion problem in implicit numerics as described e.g. in Dorfi et 
al. (2004) [ _] will remain the most time consuming part even or especially in multiple 
dimensions. 
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5 Conclusions and Prospective 



We present versatile conceptual and methodical fundamentals for strong conservation 
numerics in general non-steady curvilinear coordinates in 2D and 3D. The first chapters 
1, 2 and 3 emphasize on mathematical rigorousness and consistency which led to the 
important finding (3.6f) where we reformulated the artificial viscosity for curvilinear 
coordinates. The relevance of this finding does not come clear until general curvilinear 
coordinates are used that is why previous calculations were not necessarily affected by 
this inconsistency. 

As lateral result we suggest a dynamic equation for the self-gravitation in multiple 
dimensions by considering the specific gravitational force and the continuity equation 
(2.32). We avoid pursuing the Laplace equation in 2D and 3D where its unique nature as 
elliptical PDE in our otherwise hyperbolic problem evokes difficulties in posing adequate 
boundary conditions. Our non-static gravitation equation will have to be studied in 
future calculations whereat its imaginable applications range from stellar physics to 
galactic and cosmological computations where self-gravitation plays a role. 

We derive the equations of RHD in strong conservation form an exemplarily sketch 
them in polar and spherical coordinates in chapter 3.3. They might be understood as 
intermediate step to more applicable implementations in general curvilinear coordinates 
and play a role as reference frame a after all. 

In chapter 4 we study adaptive grid generation for curvilinear coordinates and find 
that orthogonality of coordinate lines can not be maintained for multi-dimensionally 
adaptive grids that contain the reference frames of polar and spherical coordinates. 
In this sense the polar and polar spherical take an exceptional position as orthogonal 
grids in 2D and 3D. However, we propose an interesting ansatz that satisfies our main 
postulations for an astrophysically applicable grid and study possible adoptions. For 
practical implementation one has to develop means of controling the skewness of the grid 
respectively the off-diagonal elements of the metric tensor adequately. A stringent design 
of such a technique would have exceeded the extent of this thesis and will have to be 
deferred to future investigations. Moreover one will have to contemplate feasible methods 
of generating and handling the enormous equations that emerge with non-orthogonal 
coordinates considered. Clearly, one will have to pursue automatized techniques like 
those presented with Mathematica. Concrete implementations shall also indicate at what 
point in the computations the geometric terms are to be calculated ideally. 
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6 Appendix 



6.1 Arguments and Details 
6.1.1 Symmetry of the Stress Tensor 

In section 2.1 we argue that the symmetry of the stress tensor is due to the conservation 
of angular momentum. If so, the rate of change of the angular momentum of a fluid 
px X V must equal the total applied torque. With dS oriented with normal n, let t be 
the surface force across this element and the i-th component be defined by f = T'^^rij. 



J (x X 5t(pu)) + J (xxt)-dS = 

V{t) dV{t) 

J dt{pe,jkX^u'')dV+ J EijkX^T^^nidS = 



V(t) dV{t) 



J [dt{pe^jkX^u'')+Vieijkx'T''')dV = 



v{t) 

Since the volume of our conservation law is arbitrary, we can focus on the integrand. 

dt{peijkX^v!')+eijkX^T^^^i+eijkx\iT^^ = 



ei.kT^^ = (6.1) 
The left hand side of equation (6.1) is zero, if T^^ = T^^ V j, k. 

6.1.2 Rankine-Hugoniot Conditions 

We consider a one-dimensional Riemann-problem for 1.4 with discontinuous initial con- 
dition 

={';'■ (6.2) 

dr, X > 
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6 Appendix 

and look for a solution d(x, t) for the left and the right hand side of the shock 

d{x,t) = { (6.3) 

\ dr, X > Ust 

with positive velocity Us- To do so, we split the integrals in (1.9) cleverly at the discon- 
tinuity and get for the first term 

CO oo oo oo x/us 

J J ddtjdxdt ~ J J ddt'ydxdt + j j ddt'ydxdt + 

-oo C 
oo oo 

+ J J ddtjdxdt 

Us/x 

oo 

—di J ^{x,0) dx — dr J j{x,0) dx + 

-oo 

oo 



-oo -oo 

oo oo 



+ I ^{x,x /Us){dr — dj^ dx 

oo 

7(2;, 0)do(a;) da; — / j{x, x/us)(^di — dr^ dx. (6.4) 





00 



For the flux term on the right hand side we proceed analogously and rewrite the integral 
over dx with t = x/ug. 

00 00 00 / Ust 00 

f{d)d^-fdxdt = !\ I jd^f{d)dx+ fdj{d)dx\dt 



—00 \— 00 Ust 

00 



j{Ust,t) [f{dl)-f{dr))dt 

00 

^ J 7(X, X/Us) (f{di) - f{dr)) dx (6.5) 





00 







-00 



{ddtj + f{d) dx^) dxdt = — J j{x, 0) do{x) dx + 

—00 

00 

+ J j{x,x/ns) (^ ^(^^)-^(^-) _ (di-dr)^ dx (6.6) 
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6.1 Arguments and Details 



(6.6) is weak solution of the Riemann-problem, if 

fidl) - f{dr) 



Us 

or 



(^di -dr)=0 (6.7) 



^ m)-f{dr) _ 

^ di — dr 



We call Us shock velocity and (6.8) Rankine-Hugoniot condition. 



6.1.3 DifFusion Approximation 

If the radiation field does not depend on the direction i.e. it is isotropic I^AiS = J'y^iS 7^ 
<^7,diff(n), the integration over the solid angle can be executed straightforward. Let be 
the normal vector n in local spherical coordinates. 



' 7,diff 



7,diff 



2n n 



2n n 



COS (f sin 

J J I-yBr sin 1? dii)dip = J J sin i? sin (/? sin 'd 

cos {} 



dMif = 





271" TT 







2n n 



-J J I^erGr sin -i? d'ddif ~ ~ J J d-ddip ■ 



/ COS (p^ sin cos tp sin ip sin cos p cos '& sin 

cos if sin if sin t?^ sin ip^ sin cos i? sin 93 sin ■(? 
\ cos (/? cos ?? sin ?? cos 'd sin (/? sin ■!? cos f?^ 



11 
c 



/ 47r 
' 3 



V 





3 ^ 

— 



1 



(6.9) 



6.1.4 Strong Conservation Form on Adaptive Grids 

In section 3.2.3 we outlined the strong conservation form for time dependent coordinates. 
Following Thompson, Warsi, Mastin [\'.\] we present the full deduction of this conclusion. 

The time derivative of a physical quantity d in the non-steady coordinate system S(^) 
relative to a coordinate system is given by i9td(^) = i9td(Q,) + V(a)d(9tX(^) which is 
basically the Reynolds transport theorem. In system (a) the time derivative of d yields 
9td(Q,) = 5fd(^) — xV(Q,)d. With x we introduced the grid velocity, the second term 
on the right hand side we also call grid advection. An advection term in a hyperbolic 
conservation law as we consider them usually gains the form 

K = dt + d\Y {ud) (6.10) 
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respectively in the transformed system 

= - X • Vd + div (ud). 



(6.11) 



Applying the strong conservation form of our spatial differential operators we developed 
in section 3.2.2, this convection term yields 



= d. - X • 



1 



A 



|g|eM 



+ 



1 



|g|u • eM 



(6.12) 



Since we deal with non-steady coordinates both the base vectors and the metric are 
functions of time. We resort to an identity from tensor analysis for the metric tensor to 
consider the temporal derivative of \/|g[ 



d. 



dt 



ei • (e2 X es) 



dtei ■ (e2 X es) + dte2 • (es x ei) + . . . 
|g|<9tei -e* 



(6.13) 



usmg e 
vector 



-1/2, 



Cj X ej.). Moreover, we remember the definition of the covariant base 



&. 



X 



dx 



(6.14) 



(a) 



which leads to an interesting conclusion for non-steady coordinates respectively adaptive 
grids. 

dt\l\^\ = \Agl^i^ • 6* (6.15) 

The variation of the volume element in time hence is connected with the derivatives of 
the grid velocity in all directions. Again we mention that even for solely radially adaptive 
grids the term x contains angular contributions, to wit x 7^ r. 

We consider following expression in order to asses the temporal change of a geometri- 
cally weighted in comparison to a plain density function d 



1 



-A 



|g| 



|g|d 



|g| 



^i^^lgld + y |g|5td 
g\di± ■ eM + \ \g\dtd 



and obtain 



dtd 



dt \/ |g|d — 9,;x • eM. 



(6.16) 



We compare that result with our convection term (6.12) to gain the strong conservation 
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form of the advection term in non-steady coordinates. 



vlg| 

= 



|g|d 



|g|i^ = dt J\g\d 



|g|d 



v|g| 

1 



-A 



|g| 



|g|x • eM 



|g|eM 
1 



|g|u • eM 



|g|u • eM 



:|e* • (u - x)d 

If we define the contravariant velocity components relative to the moving grid by 

U' = e' • (u - x) 

the upper equations yields 



(6.17) 



\g\K = dt J\g\d +di J\g\U'd 



(6.18) 
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6.2 Mathematica Files 

On the following pages some examples of Mathematica Notebooks are presented that were 
used for tensor analytical computations in several coordinates and generation of strong 
conservation form equations. The analysis of the grid generation toy models presented 
in sections 4.4 and 4.6 conclude the Appendix. 

Technical remark: in order to ensure that integral respectively conserved quantities 
are not divided automatically by the computer algebra system, we make use of Patterns. 
In this sense the partial derivatives act on braced expressions like 9j [. . . ] that would be 
processed for the discretization. 

Spherical Coordinates (spherical_coordinates .nb) Differential geometric definitions 
and relations for base vectors, metric components and Christoffel symbols, exem- 
plarily executed with spherical coordinates. If the initial definition of the coordi- 
nate transformation gets modified to a nonorthogonal coordinate system, minor 
changes have to be made. 

Artificial Viscosity in Spherical Coordinates (artificial_viscosity.nb) Using results 
from spherical_coordinates .nb we compute co- and contravariant components 
of the viscous pressure tensor in spherical coordinates. Symmetry of Q and van- 
ishing trace are evaluated explicitly. 

Equation of Motion in Spherical Coordiantes (E0M_3D_sph.nb) Using the previous def- 
initions for geometric quantities and the viscous pressure we generate the equation 
of motion in spherical coordinates in strong conservation form for non-steady grids. 

Polynomial Coordinates Ansatz (poly_coordinates .nb) Graphical analysis of a poly- 
nomial ansatz for orthogonal family of curves. 

Polar Non-Orthogonal Grid Ansatz (quasi_polar_grid.nb) Generation of a non-orthogonal 
grid that in principle satisfies all postulations from 4.2 besides orthogonality. With 
the coordinate transformation determined one could go run through upper files and 
generate strong conservation form equations on non-steady non-orthogonal grids. 
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Spherical Coordinates 

Differential Geometric Definitions and Relations 



■ Coordinate Transformation 

ln[1]:= X := rSin[theta] Cos [phi] 
y := rSin[theta] Sin [phi] 
z : = r Cos [theta] 

ln[4]:= lij = {{D[x, r] , D[x, theta], D[x, phi]}, {D[y, r] , D[y, theta], D[y, phi]}, 
{D[z, r] , D[z, theta], D[z, phi]}}; (* Jacobian of transformation *) 

ln[5]:^ li jtransposed : = Transpose[lij] 

■ Covariant Components of Metric Tensor 

ln[6]:= gij = li jtransposed. lij // FullSlmplify (* from gij=ei.ej *) 

Out[6]= {{1, 0, 0}, {O, O}, {O, 0, Sin [theta] ^ } } 

ln[7]:= g = Det[gij] 
Out[7]= sin [theta] ^ 

ln[8]:^ rootg = Assi]ming[ {r > 0, Sin[theta] >0}, Ref ine [Sqrt [g] ] ] 

Out[8]= Sin [theta] 

■ Covariant Base Vectors 

ln[9]:^ er = 11 jtransposed[ [1] ] (* tt-st row of Jacobian defines the tt-st convariant base vector *) 

Out[9]= {Cos [phi] Sin [theta]. Sin [phi] Sin [theta], Cos [theta] } 

ln[10]:= etheta = 11 jtransposed [ [2] ] 

Out[10]= {rCos[phi] Cos [theta] , rCos[theta] Sin[phi], - r Sin [theta] } 
ln[11]:= ephi = lijtransposed[ [3] ] 

Out[11]= {-rSin[phi] Sin[theta], r Cos [phi] Sin[theta], 0} 

■ Contravariant Base Vectors 

ln[12]:= Er = 1 / rootg * Cross [etheti, Sphi] // FullSlmplify 

(* this identity is only satisfied for orthogonal coordinates; 
for nonorthogonal coordinats one needs to raise indices g^'^ej *) 

Out[12l= {Cos [phi] Sin [theta]. Sin [phi] Sin [theta], Cos [theta]} 



2 I spherical_coordinates.nb 



ln[13]:= Etheta = 1 / rootg * Cross [Cphi , e,] // FullSimplify 

r Cos [phi] Cos [theta] Cos [theta] Sin [phi] Sin [theta] 

Out[13]= \ , , 

^ r r r 

In[l4]:= Ephi = 1 / rootg * Cross [e^, etheta] // FullSimplify 

Csc [theta] Sin [phi] Cos [phi] Csc [theta] 



Out[14]= 



■ Covariant, Mixed and Contravariant Components of the Metric Tensor 

ln[15]:= Table [inetricg[i, j] =ei.ej //FullSimplify, {i, {r, theta, phi}}, {j, {r, theta, phi}}] 

Out[15]= {{1, 0, 0}, {O, O}, {O, 0, Sin [theta] ^ } } 

ln[16]:= Table [Metricg[i, j] =Ei.Ej // FullSinplify, {i, {r, theta, phi}}, {j, {r, theta, phi}}] 

f r 1 1 r Csc [theta] ^ -, -, 

0U.I16,. {{1, 0, 0}, {O, -, O}, {O, 0, }} 

ln[17]:= Table [delta [i, j] =Ei.e] // FullSiit5)lify, {i, {r, theta, phi}}, {j, {r, theta, phi}}] 

Out[17]= {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}} 

■ Christoffel Symbols 

ln[l8]:^ Table [Table [Table [Gammer [sigma, mu, nu] = 

1 / 2 * Sum [ (Ejigma ■ Erho) * (D [ (e„u . erho) f mu] + D [ (e„u . erho) / nu] - D [ (e„u . e„u) , rho] ) , 
{rho, {r, theta, phi}}] // FullSiii5)lify, 
{sigma, {r, theta, phi}}], {mu, {r, theta, phi}}], {nu, {r, theta, phi}}] 

Out[18]= {{{0, 0, 0}, {o, -, o|, {o, 0, -||, {{o, -, o|, {-r, 0, 0}, {0, 0, Cot [theta] }} , 
1 

||o, 0, — J, {0, 0, Cot[theta]}, {- r Sin [theta] ^ , -Cos[theta] Sin[theta], 0}|| 



Artificial Viscosity 

Viscous Pressure Tensor 



■ Definitions of Viscous Pressure 

In order to proof consistency of our modified artificial viscosity we use different definitions and relations. 

In[19]:= Table [ 

qschlange[i, j] = 

1/2* (d[i] [u[j]] +d[j] [u[i]] - Sum [ Gammer [ k , i, j] *u[k], {k, {r, theta, phi}}] - 

Sum[Gammer [k, j, i]*u[k], {k, {r, theta, phi} } ] ) - 1 / 3 * metricg[i, j]*div[u], 
{i, {r, theta, phi}}, {j, {r, theta, phi}}]; (* covariant con^onentes *) 

ln[20]:= Table [qschlangeoben [m, n] = Sirni [ (Metricg [m, k] *Metricg[n, 1]) * 

(1 / 2 * (d[k] [u[l] ] + d[l] [u[k] ] - Sum [ Gammer [ o , k, 1] *u[o], {o, {r, theta, phi}}] - 
Sum[Gammer [o, 1, k] *u[o], {o, {r, theta, phi}}]) - 
1 / 3 * metricg[m, n] *div[u]), {k, {r, theta, phi}}, {1, {r, theta, phi}}], 
{m, {r, theta, phi}}, {n, {r, theta, phi}}]; (* contravariant components *) 

ln[21]:= Table [ 

qoben[m, n] = Sum[ (Metriog[m, i] *Metriog[n, j]) qschlange[i, j] , 

{i, {r, theta, phi}}, {j, {r, theta, phi}}], 
{m, {r, theta, phi}}, {n, {r, theta, phi}}]; (* contravariant con^onents *) 

ln[22]:= Table [ 

qinix[i, j] = Sum [Metricg [ 1 , i] * qschlange [1, j], {1, {r, theta, phi}}], 
{i, {r, theta, phi}}, {j, {r, theta, phi}}]; (* mixed con^onents *) 

ln[23]:= Table [ 

qinix2[i, n] = Sum[metricg[ j, n] * qschlangeoben [ i , j], {j, {r, theta, phi}}], 
{i, {r, theta, phi}}, {n, {r, theta, phi}}]; (* mixed conponents *) 
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■ Components Explicitly 

ln[24]:= Table [qmix[m, n] , {m, {r, theta, phi}}, {n, {r, theta, phi}}] // FullSimplify // Expand 

«div[u] u [theta] 1 1 
+ d[r] [u[r] ] , + -d[r] [u [theta] ] + -d [theta] [u[r] ] , 
3 r 2 2 

u[phi] 1 1 

+ -d[phi] [u[r]] + -d[r] [u[phi]] , 

r 2 2 J 

|- u[theta] d[r] [u [theta] ] d [theta] [u [r] ] div[u] u[r] d [theta] [u [theta] ] 
L J.3 2 2 3 r 

Cot [theta] u [phi] d [phi] [u [theta] ] d [theta] [u [phi] ] 



}' 



r-" 2 2 

Csc[theta]^ u[phi] Csc[theta]^ d[phi] [u[r] ] Csc[theta]^ d[r] [u[phi] ] 



2 2 

Cot [theta] Csc [theta] ^ u [phi] 



Csc[theta]^ d[phi] [u[theta] ] Csc[theta]^ d[theta] [u[phi] ] 



2 2 
div[u] u[r] Cot [theta] u [theta] Csc [theta] ^ d [phi] [u [phi ] ] 



}} 



ln[25]:= Do [Print [Q[m, n] , " = ", qschlange [m, n] // FullSinplify // Expand] , 
{m, {r, theta, phi}}, {n, {r, theta, phi}}] 

div [u] 

Q[r, r] = + d[r][u[r]] 

3 

u [theta] 1 1 

Q[r, theta] = + - d[r] [u [theta] ] + — d [theta] [u[r] ] 

r 2 2 

u[phi] 1 1 

Q[r, phi] = + -d[phi] [u[r]] + -d[r] [u[phi]] 

r 2 2 

u [theta] 1 1 

Q [theta, r] = + - d[r] [u [theta] ] + — d [theta] [u[r] ] 

r 2 2 

1 

Q [theta, theta] = div [u] + r u [ r] + d [theta] [u [theta] ] 

3 

1 1 

Q[theta, phi] = -Cot[theta] u[phi] + — d[phi] [u[theta] ] + — d[theta] [u[phi] ] 

2 2 

u[phi] 1 1 

Q[phi, r] = + -d[phi] [u[r]] + -d[r] [u[phi]] 

r 2 2 

1 1 

Q [phi, theta] = -Cot [theta] u [phi] + — d[phi] [u [theta] ] + — d [theta] [u [phi] ] 

2 2 

Q[phi, phi] = 
1 

div [u] Sin [theta] ^ + r Sin [theta] ^ u [r] + Cos [theta] Sin [theta] u [theta] + d[phi] [u [phi] ] 

3 

ln[26]:= Do [Print [Qcontra[m, n] , " = ", qschlangeoben [m, n] // FullSimplify // Expand] , 
{m, {r, theta, phi}}, {n, {r, theta, phi}}] 
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div [u] 

Qcontra[r, r] = + d[r][u[r]] 

3 

u[theta] d[r] [u [theta] ] d [theta] [u [r] ] 
Qcontra[r, theta] = + — 



Qcontra[r, phi] = 
Qcontra [theta, r] 
Qcontra [theta, theta] 



2 r'' 2 
Csc[theta]^ u[phi] Csc[theta]^ d[phi] [u[r] ] Csc[theta]^ d[r] [u[phi] ] 



2 2 
u[theta] d[r] [u [theta] ] d [theta] [u [r] ] 



J.3 2 2 

div[u] u[r] d[theta] [u[theta] ] 



3 r2 



Qcontra [ theta, phi] = 

Cot [theta] Csc [theta] ^ u [phi] Csc [theta] ^ d[phi] [u [theta] ] Csc [theta] ^ d [theta] [u [phi] ] 



2 r" 2 r" 

Csc[theta]^ u[phi] Csc[theta]^ d[phi] [u[r] ] Csc[theta]^ d[r] [u[phi] ] 



Qcontra[phi, r] 

r-^ 2 2 

Qcontra [phi, theta] = 

Cot [theta] Csc [theta] ^ u [phi] Csc [theta] ^ d[phi] [u [theta] ] Csc [theta] ^ d [theta] [u [phi] ] 



2 r" 2 r" 

Csc[theta]^ div[u] Csc[theta]^ u[r] 



Qcontra [phi, phi] 

3 r'^ 

Cot[theta] Csc[theta]^ u[theta] Csc[theta]'' d[phi] [u[phi] ] 



■ Symmetry 

In fact Q is symmetric in its co- and contravariant components. 

In[27]:^ Table [qschlangeoben[m, n] = = = qschlangeoben[n, m] , {m, {r, theta, phi}}, {n, {r, theta, phi}}] 

Out[27]= {{True, True, True}, {True, True, True}, {True, True, True}} 

ln[28]:= Table [qschlange [m, n] === qschlange [n, m] , {m, {r, theta, phi}}, {n, {r, theta, phi}}] 

Out[28]= {{True, True, True}, {True, True, True}, {True, True, True}} 

ln[29]:= Table [qinix[m, n] === qniix[n, m] , {m, {r, theta, phi}}, {n, {r, theta, phi}}] 

Out[29]= {{True, False, False}, {False, True, False}, {False, False, True}} 



■ Trace of Q 



The trace vanishes as expected. 
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ln[30]:= Sum[qinix[m, m] , {m, {r, theta, phi}}] 

div[u] 1 , f 1 , , 
Out[30]= + — Csc[theta] r div [u] Sin [theta] + 

3 r2 [3 

— (2 r Sin [theta] ^ u [r] + 2 Cos [theta] Sin [theta] u [theta] + 2 d[phi] [u[phi] ] ) J + 



d[r] [u[r] ] + 



J div[u] + ^ (2 r u[r] +2 d[theta] [u[theta] ] ) 



ln[31]:= Sum[<ynix[m, m] , {m, {r, theta, phi}}] / 
{div[u] 



2u[r] Cot [theta] u [theta] Csc [theta] ^ d[phi] [u [phi] ] 



r r" 



d[theta] [u [theta] ] 
d[r] [u[r]] + } // FullSiii?>lify // Expand 



Out[31]= 



ln[32]:= Sum [qmix2 [m, m] , {m, {r, theta, phi}}] /. 

2u[r] Cot [theta] u [theta] Csc[theta]^ d[phi] [u[phi] ] 



{div[u] 



r 



d[theta] [u [theta] ] 
d[r] [u[r]] + } // FullSiii?>lify // Expand 

__5 ' 



Out[32]= 



Equation of Motion in Spherical Coordinates 

Generation of Strong Conservation Form Equation 



■ Definition of Contravariant Components of Tensors 

Table [rschlange [1, j] = rho * (u[i] -xdot[i]) *u[j] // FullSimplify, 

{i, {r, theta, phi}}, {j, {r, theta, phi}}] 
Table[p[i, j] = p * Metricg[i, j] //FullSimplify, {i, {r, theta, phi}}, {j, {r, theta, phi}}] 
Table[q[i, j] = Siaiii[ (Metricg[i, m] *Metricg[j, n] ) 

(1 / 2 * (d[n] [u[m] ] + d[n] [u[m] ] - Sum[Gammer [k, m, m] * u[k] , {k, {r, theta, phi} } ] - 

Svim [Gammer [k, m, n] *u[k], {k, {r, theta, phi} } ] ) - 1 / 3 * metricg[n, m] * div[u] ) , 
{n, {r, theta, phi}}, {m, {r, theta, phi}}] // FullSiii?>lify, 
{i, {r, theta, phi}}, {j, {r, theta, phi}}] 
(* note that there emerge covariant velocity conponents in Q due to its definition *) 

Out[33]= {{rhou[r] (u[r] - xdot [r] ) , rho u [theta] (u[r] -xdot[r]), rho u [phi] (u[r] -xdot[r])}, 
{rhou[r] (u [theta] - xdot [theta] ) , rho u [theta] (u [theta] - xdot [theta] ) , 
rhou[phi] (u [theta] - xdot [theta] )} , {rhou[r] (u[phi] -xdot [phi]), 
rho u [theta] (u[phi] - xdot [phi]), rho u [phi] (u[phi] - xdot [phi] )} } 



i- (- p , (- pCsc [theta] -, -, 
Out[34]= |{p, 0, 0}, jo, — , Oj, jo, 0, ^ II 



|- |- div[u] 
Out[35]= II + d[r][u[r]]. 



u [theta] - 2 r d [theta] [u [r] ] Csc [theta] ^ (u[phi] - 2 r d[phi] [u[r] ] ) ^ 



2 r^ 2 r^ 

,2 



■ r^ u[r] - u [theta] + 2 r d[r] [u [theta] ] -r^ div[u] +3 (r u[r] + d[theta] [u [theta] ] ) 



2 r^ 3 r 

Csc [theta] ^ (-Cot [theta] u [phi] + r u [r] +2 d[phi] [u [theta] ] ) 



2 r" 



}' 



. r^ u[r] + r Cot [theta] u [theta] - Csc [theta] ^ (u [phi] - 2 r d[r] [u[phi] ] ) 

2 

-Cot [theta] Csc [theta] ^ u [phi] + r u [r] + Cot [theta] u [theta] + 2 Csc [theta] ^ d [theta] [u [phi] ] 

2 r" 

Csc [theta] ^ (-r^ div[u] +3 (r u[r] + Cot [theta] u [theta] + Csc [theta] ^ d[phi] [u [phi] ] ) ) 



1) 



■ Equation of Motion in Spherical Coordinates and Strong Conservation Form 



ln[36]:= Table [d[t] [Sum[Sqrt [g] *rho*u[n] *e„[[k]], {n, {r, theta, phi}}]] +Sum[ 

Sum[d[i] [Sqrt [g] *( rschlange [i, j]+p[i, j]+Qkomp[i, j] ) * ej [ [k] ] ] , {i, {r, theta, phi} }] , 
{j, {r, theta, phi}}] + rho * Sum[d[l] [Sqrt [g] * Phi * El [ [k] ]] , {1, {r, theta, phi}}] - 
4 Pi /c* Sqrt [g] * {s™[H[o] *eo[[k]], {o, {r, theta, phi}}]) /. 
{Sin[theta] ^ Sqrt [1 - mu^2] , d[theta] ^ - Sqrt [ 1 - mu" 2] d[mu] , 
Csc[theta] -» 1 / Sqrt [1 - mu^ 2] , Cos[theta] ^ mu} , {k, {1, 2, 3}}] 
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Out[36]= 



4/7-^(1 -mu^) - mu^ Cos [phi] H[r] + mu r Cos [phi] H[theta] - ^1 - mu^ rH[phi] Sin [phi] 

c 

|- ^1 - mu^ d[mu] j [-^1 - mu^ r ^(l -mu^) r" Sin[phi] 

(Qkomp [theta, phi] +rhou[phi] (u[theta] - xdot [theta] ) ) ] + a^I - mu^ d[mu]j[ 
^1 - mu^ ^(l - mu^) Cos [phi] ( Qkomp [theta, r]+rhou[r] (u [theta] - xdot [theta] ))] + 



— + Qkomp [theta, theta] +rhou [theta] (u [theta] - xdot [theta] ) j j +d[phi] 



/1-mu^ r 



1-mu^) Sin [phi] + Qkomp [phi, phi] +rhou[phi] (u[phi] - xdot [phi]) + 

1^(1 -mu2) r2 jJ 

d[phi] l^^l -mu^ ^ ( 1 - mu^ ) Cos [phi] (Qkomp[phi, r]+rhou[r] (u [phi] - xdot [phi] )) j + 

d [phi ] l^mu r ( 1 - mu^) r" Cos[phi] (Qkomp[phi, theta] + rho u [theta] (u[phi] - xdot [phi] )) j + 



rho 



1 - mu d [mu] 



_ mu Phi ^(l - mu^) r" Cos [phi] 



Phi , / (1 - mu^) Sin[phi] i , 

d[phi] ^ — ^ +d[r] U/l-mu^ PhiW(l-mu=^) r" Cos[phi] 

/ 1 - mu^ r 



d[r]^-^l-mu^ r^ll-mu^jr" Sin[phi] (Qkomp [r, phi] + rho u [phi] (u [r] - xdot [r] ) ) J + 
d[r] [^1 - mu^ -^(l -mu^) r" Cos [phi] (p + Qkomp[r, r] +rhou[r] (u[r] - xdot [r] ) ) ] + 
d[r] l^mu r ^(l -mu^) Cos[phi] (Qkomp[r, theta] +rhou[theta] (u[r] - xdot [r] ) ) j + 
d[t] [-^1 -mu^ r^(l-mu^) r" rhoSin[phi] u[phi] + 

^1 - mu^ ^ (l -mu^) r" rho Cos [phi] u [r] + mu r ^(l - mu^) rho Cos [phi] u [theta] ] , 
4 71 -^(l - mu^) r* [\Jt- - mu^ r Cos [phi] H[phi] + ^yT^^niu^" H[r] Sin [phi] +murH [theta] Sin [phi] 

c 

^- ^1 - mu^ d[mu] j - mu^ r ^ (l - mu^) r^ Cos [phi] 
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(Qkomp [theta, phi] +rhou[phi] (u[theta] - xdot [theta] ) ) j + ^-^1 - mu^ d[mu]j|^ 



1 - mu ^/ 1 - mu r Sin [phi] (Qkomp [theta, r] + rho u [r] (u [theta] - xdot [theta] )) + 



- mu^ d[mu] j Jmu r ^(l - mu^) r* Sin[phi] 



f — + Qkomp [theta, theta] + rho u [theta] (u [theta] - xdot [theta] ) ) | + d [phi] \ ^Jl - mu^ r 



^(l - mu^) r" Cos [phi] 



P 

+ Qkomp [phi, phi] + rho u [phi] (u[phi] - xdot [phi]) 

2\ ^2 



1 - mu ) r 



d [phi ] l^mu r ( 1 - mu^) r^ Sin[phi] (Qkomp[phi, theta] +rhou[theta] (u[phi] - xdot [phi] )) j + 

I \ mu Phi W(l - mu^) r" Sin [phi] 

yi-mu^ d[mu]j[ ^ J + 



rho 



Phi ./(I -mu^) r" Cos [phi] i , 

d[phi] ^ — ^ +d[r]Uyi-mu^ Phi W ( 1 - mu^ ) r" Sin[phi] 

/ 1 - mu^ r 



d [ r ] J ^Jl-nva^ r ^ ( 1 - mu^ ) Cos [phi] (Qkomp [r, phi] +rhou[phi] (u[r] - xdot [ r ] ) ) j 
d [ r ] I 1 - mu^ \J ^ } Sin[phi] (p + Qkomp [r, r] +rhou[r] (u[r] ~ xdot [ r ] ) ) j + 



d[t] [^1 - mu^ r ^(l - mu^) r" rho Cos [phi] u [phi] 



-yl-mu^ -^J i^l - mu^ j r^ rho Sin [phi] u[r] + mu r 



r rho Sin [phi] u [r] + mu r ^/ 1 - mu ) r rho Sin [phi] u [theta] , 



71 -^(l - mu^) r^ (mu H[r] - ^1 - mu^ rH[theta] 



+ -a/1 -mu-^ d[mu] [0] + 



^-^1 - mu^ d[mu] j |^mu ^(l - mu^) r^ (Qkomp [theta, r] +rhou[r] (u[theta] - xdot [theta] )) j + 

|- a/i - mu^ d[mu] j [ 

-\Jl - mu^ r ^(l - mu^) r^ ^ — + Qkomp [theta, theta] + rho u [theta] (u [theta] - xdot [theta] ) j j + 
d[phi] [0] + d[phi] |^mu ^ ( 1 - mu^ ) r^ (Qkomp [phi, r] + rho u[r] (u[phi] - xdot [phi] ) ) ] + d[phi] [ 
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rho 



^1 - mu^ r^(l-mu^) r" (Qkomp[phi, theta] +rhou[theta] (u [phi] - xdot [phi] ) ) ] + d [ r ] [ ] + 

1 - mu^ Phi ^(l -mu^) r 



|- ^1 - mu^ d[mu] J 



+ d[phi] [0] + d[r] [mu Phi ^(l -mu^) r" 



d [ r ] |mu ^ 1 1 - mu^ J (p + Qkomp [r, r] + rhou[r] (u[r] - xdot[r]))j 



d [ r ] [ " ^ 1 ~ 't'U^ r ^ ( 1 - mu^) (Qkomp [ r , theta] + rho u [theta] (u [r ] - xdot [r ] ) ) J + 
d[t] l^mu ^ ( 1 - mu^ ) r" rhou[r] - ^1 - mu^ r^(l-mu^)r'' rho u [theta] ] J 



■ EOM with Q explicitly 

ln[37]:= Table [d[t] [Siaiii[Sqrt [g] *rho*u[n] *e„[[k]], {n, {r, theta, phi}}]] + 

Sum[Sum[d[i] [Sqrt[g] * ( rschlange [ i , j] +p[i, j] +q[i, j]) *e][[k]]], {i, {r, theta, phi}}], 
{j, {r, theta, phi}}] - +rho * Sum[d[l] [Sqrt [g] * Phi * Ei [ [k] ] ] , {1, {r, theta, phi}}] 
4 Pi / c * Sqrt[g] * (Sum[H[o] *eo[[k]], {o, {r, theta, phi}}]) /. 
{Sin[theta] -» Sqrt [1 - mu'' 2] , d[theta] -» - Sqrt [1 - mu* 2] d[mu] , 
Csc[theta] ^ 1 / Sqrt [1 - mu'^2] , Cos [theta] -» mu} , {k, {1, 2, 3}}] 



Out[37]= 



l-^l - mu^ d[mu] j [mu r ^(l -mu^) r^ Cos [phi] 



■ + rho u [theta] (u [theta] - xdot [theta] ) + 



r^ div[u] +3 r u[r] + - mu^ d[mu] [u[theta] ] 



3 r« 



l^-^Jl- mu^ d[mu] j [-^1 - mu^ r ^(l -mu^) r" Sin[phi] 



rhou[phi] (u[theta] - xdot [theta] ) + 



-Cot [theta] u[phi] + r u[r] +2 d[phi] [u [theta] ] 



|- ^1 - mu^ d[mu] j [^1 - mu=^ ^ ( 1 - mu^ ) 



rhou[r] (u [theta] - xdot [theta] ) 



2 (l -mu^) 
r" Cos [phi] 

r^ u[r] - u [theta] + 2 r d[r] [u [theta] ] 



2 r^ 
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d[phi] l^mu r ^(l - mu^) r" Cos [phi] 



rhou[theta] (u[phi] -xdot[phi]) + 



Cot [theta] u[phi] 



+ ru[r] + Cot [theta] u [theta] + 



[4 



2 -^/l-mu^ d[mu] [u[phi] 



2 r" 



d[phi] - mu^ r -^(l -mu^) r" Sin[phi] 



1 1 - mu^ J 



■ + rhou[phi] (u[phi] -xdot[phi]) + 



-r^ div[u] +3 ru[r] + Cot [theta] u[theta] + 



d [phi ] [u [phi] 



3 1 - mu^ r" 



d[phi] [^1 - mu^ ^(l - mu^) Cos [ph 



i] rhou[r] (u[phi] -xdot[phi]) + 



u [r] + r Cot [theta] u [theta] 



u [phi] - 2rd[r][u[phi]] ~ 
1-mu^ 



2 



] 4 TT •^y ( 1 - mu^ ) 



rho 1^1 -mu^ Cos [phi] H[r] + mu r Cos [phi] H [theta] -^l-mu"^ rH[phi] Sin [phi] 
mu Phi ^( 1 - mu^) Cos [phi] 



|- ^1 - mu^ d[mu]jj- 



■ d[phi] 



Phi ^/ 1 - mu'' Sin[phi] 



■I 



1 - mu r 



rho u [theta] (u[r] - xdot [r] ) 



u [theta] - 2r -A/l-mu d[mu] [u[r] 



2 



- ^Jl-rny^ r ^ ( 1 - mu^ ) 



1 - mu r A / [ 1 - mu ) r Sin[phi] rhou[phi] (u[r] - xdot [r] ) 



-d[r] 



u[phi] - 2 r d[phi] [u[r] ] 



2 1 - mu^ 
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d[r] 



[a/i - mu^ ^(l -mu^) 



div[u] 

r" Cos[phi] |p + rhou[r] (u [ r ] - xdot [r ] ) + d [r ] [u [ r ] ] + 



d[t] I - ^1 - mu^ T ^ 1^1 - mu^ J rhoSin[phi] u[phi] + 
-X^l - mu^ ^( 1 - mu^) rho Cos [phi] u[ 



1 - mu d[mu] mu r a / 1 - mu r Sin [phi] 



■ + rho u [theta] (u[theta] - xdot [theta] ) + 



div[u] +3 r u[r] + -a/1 ' d[mu] [u[theta] ] 



3 r« 



|- a/i - mu^ d[mu] j l^^/r^^mu^ r ^(l - mu^) Cos [phi] 



rho u [phi] (u [theta] - xdot [theta] ) + 



i^-^Jl ~ mu^ d[mu] j [^/l " mu^ -^(l - mu^ ) r' Sin [phi] 



-Cot [theta] u [phi] + r u [r] +2 d[phi] [u [theta] ] 



2 (l - mu^) r* 



rhou[r] (u [theta] - xdot [theta] ) + 



r^ u [r] - u [theta] + 2 r d [r] [u [theta] ] 



2 r^ 



d[phi] mu r / 1 - mu ) r Sin[phi] 



rho u [theta] (u[phi] - xdot [phi]) + 



Cot [theta] u[phi] 



+ ru[r] + Cot [theta] u [theta] + 



2 \—Jl-inu' d[muj [u[phl] 



2 r" 



d[phi] - mu^ r ^(l - mu^) r^ Cos [phi] 



■ + rho u [phi ] (u[phi] - xdot [phi]) + 



-r^div[u] +3 ru[r] + Cot [theta] u[theta] + 



3 1 - mu" r' 



1 - mu j r 

d[phi] [u[phi] ] \ \ 
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d [phi - mu^ ( 1 - mu^ j Sin[phi] rhou[r] (u [phi] - xdot [phi] ) 



u [r] + r Cot [theta] u [theta] 



u[phi]-2 r d[r] [u[phi] ] '\ 
l-mu^ 



2 



] - -4 TT ( 1 - mu^ ) 



rho 1 - mu r Cos [phi] H[phi] + ^ 1 - mu H[r] Sin[phi] + mu r H [theta] Sin[phi] 

''/ I \ mu Phi , /(l - mu^) r" Sin [phi] Phi , / ( 1 - mu^) Cos [phi] 

l^-^l-mu^ d[mu]j[ ^ J+d[phi][ ^ 



4 



1 - mu r 



rho u [theta] (u[r] - xdot [r] ) 



u [theta] - 2 r - A/l - mu d [mu] [u [ r ] ] 



2 



-d[r] 



l-mu^ rW(l-mu^)r* Cos [phi] 



rho u [phi] (u[r] -xdot[r]) 



u[phi] - 2 r d[phi] [u[r] ] 



2 1 - mu=^ r^ 



d[r] [^1 - mu^ ^(l -mu^) r" Sin[phi] p ^ + rhou[r] (u [ r ] - xdot [r ] ) + d [r ] [u [ r ] ] 



d[t] ^-^1 - mu^ r^(l-mu^) r^ rho Cos [phi] u[phi] +-\jl-mu^ ^^(l - mu^ ) 



r rho Sin [phi] u[r] + 



i^-'\Jl - mu^ d[mu] j [^-^1 - mu^ r ^(l - mu^) r* 



■ + rho u [theta] (u [theta] - xdot [theta] ) + 



r^ div[u] +3 [r u[r] + f- ^1 - mu^ d[mu] 1 [u[theta] ] 1 



3 r« 



] + ^-^1 -mu2 d[mu] j [ 



mu -I / 1 - mu ) r rhou[r] (u[theta] - xdot [theta] ) + 



r^ u[r] - u [theta] + 2 r d[r] [u [theta] ; 



2 r^ 
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d[phi] [0] +d[phi] [-^1 -mu= 



l-mu=^) r" 



rhou[theta] (u[phi] -xdot[phi]) + 



Cot [theta] u[phi] 



+ ru[r] + Cot [theta] u [theta] + 



[4 



2 -^/l-mu^ d[mu] [u[phi] 



2 r" 



d[phi] [mu ^(l - mu^) 



rhou[r] (u[phi] - xdot [phi] ) + 



u [r] + r Cot [theta] u [theta] 



u[phi]-2 rd[rj [u[phij ] 



2 



^) rho |mu H[r] - ^ 



\ Jl-mu^ Phi J(l -mu2) r" ^ 

1 -mu=^ d[mu] I [ J +d[phi] [0] +d[r] |^mu Phi y (l -mu^) r" 



d[r] [-a/i -mu^ r ^(l -mu^) r" 



rho u [theta] (u[r] - xdot [r] ) 



u[theta] -2 r -a/I " ""U^ d[mu] [u[r] ] 



2 



d[r] [mu sj [l -rau^) r" (p - ^^^^"^ + rho u[r] (u[r] - xdot [r] ) + d[r] [u[r]] ]] + 



d[t] [mu a/ (1 - mu^) r' rho u [r] - ^/l - mu^ r ^ ( 1 - mu^) r" rhou[theta] 



Polynomial Coordinates Ansatz 

Grid Generation via Orthogonal Curves 



■ Ansatz for Polynomial Curves 

ln[1]:= I7[x_, r_, t_] = aO [r, t] + al [r, t] * x + a2 [r, t] * x''2 + a3 [r, t] * x'-B 
(*+a4 [r, t] *x*4+a5 [r,t] *x'>5*) ; 

First derivative to determine slope 

ln[2]:= g[x_, r_, t_] = D[r7[x, r, t] , x] ; 

Second derivative 

ln[3]:= h[x_, r_, t_] =D[r;[x, r, t] , x, x] ; 



■ Rotation lUlatrix to Rotate Coordinate System 

Since infinite slope cannot be requiremeent for linear equation system, the whole coorindate system gets rotated by 45 degree. So 
the slopes are + 1 and - 1. 

In[4]:= A[a_] = {{Cos[a], -Sin[a]}, {Sin[a], Cos[a]}}; 

ln[5]:^ coef £ = Table [ 

PI = A[Pi / 4] .{0.1 + n* 0.13, 0}; 
P2 = A[Pi / 4] .{0, 0.1 + n*0.1}; 

Solve[{J7[Pl[[l]], r, t] ==P1[[2]], r7[P2[[l]], r, t] ==P2[[2]], 

g[Pl[[l]], r, t] == -1, g[P2[[l]], r, t] == 1 (*, f [0 , r , t] ==0. l+n*0 . 12, g[0, r, t] ==0*) } , 
{aO[r, t], al[r, t] , a2[r, t] , a3[r, t] (*, a4 [r, t] , a5 [r , t] *) } ] , {n, 0, 10}] 
(* generation of polynomials Intersecting with two equldistanced points at x- 
and y-axis *) 

Out[5]= {{{aO[r, t] ^ 0.106066, al[r, t] ^0., a2[r, t] ^-7.07107, a3[r, t] ^0.}}, 

{{aO[r, t] ^0.226564, al[r, t] ^0.173909, a2[r, t] ^-3.24084, a3[r, t] ^-1.5093}}, 
{{aO[r, t] ^0.346169, al[r, t] ^0.226146, a2[r, t] ^-2.08962, a3[r, t] ^-0.834794}}, 
{{aO[r, t] ^ 0.465575, al[r, t] ^ 0.251258, a2[r, t] 54026, a3[r, t] ^- 0.510661}}, 

{{aO[r, t] ^0.584904, al[r, t] ^0.266012, a2[r, t] ^-1.21921, a3[r, t] ^-0.341655}}, 
{{aO[r, t] ^ 0.704197, al[r, t] ^ 0.27572, a2 [r, t] ^- 1.00877, a3[r, t] ^- 0.243865}}, 
{{aO[r, t] ^0.823468, al[r, t] ^0.282592, a2[r, t] ^-0.860221, a3 [r, t] ^-0.182541}}, 
{{aO[r, t] ^ 0.942727, al[r, t] ^ 0.287713, a2[r, t] ->- 0. 749781, a3 [r, t] . 141659} } , 
{{aO[r, t] ^ 1.06198, al [r, t] ^ 0.291675, a2 [r, t] ~>-0. 664457, a3[r, t] ^- 0.113079}}, 
{{aO[r, t] ^1.18122, al [r, t] ^0.294833, a2 [r, t] ^-0.59656, a3[r, t] ^-0.0923307}}, 
{{aO[r, t] ^ 1.30046, al[r, t] ^ 0.297408, a2 [r, t] ^-0.541248, a3 [r, t] ^- 0.0768}}} 
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ln[6]:= Rotate [ 

Show [Plot [ {Evaluate [7; [x, r, t] /. {coef f} ](*, Evaluate [f [x, r,t] /. {ergmore2} ]*) , Abs [x] , 0}, 
{x, -1.2, 1.2}, PlotRange -» {-0.1, 1.5}, AspectRatio ^ Automatic, Background ^ White, 
Axes -> None, PlotStyle -» {{Black, Thick}, {Gray, Thick}, {Gray, Thick}}] (*, 
Graphics [{Thick,Red,Arrow[{ {0,0}, {1, 0} }], Thick, Blue, Arrow [{{ 0, 0} , {0,1}}], 
Thick, Black, Dashed,Arrow[{{0,0},l*{l,0}+l*{0,l}}]}]*) ] , -45 Degree] 



Out[6]= 




Orthogonal Curves 



ln[71:= 



ln[15]:= 



pol[x_] 


= f)[x. 


r, t] /. {coeff[[l]]} 


(* define polynomial 


po2 [x_] 


= ^[x. 


r, t] /. {coeff[[2]]} 




po3 [x_] 


= n[x. 


r, t] /. {coeff[[3]]} 




po4 [x_] 


= n[x. 


r, t] /. {coeff[[4]]} 




po5 [x_] 


= f)[x. 


r, t] /. {coeff[[5]]} 




po6 [x_] 


= f)[x. 


r, t] /. {coeff[[6]]} 




po7[x_] 


= ^[x 


r, t] /. {coeff[[7]]}; 


po8 [x_] 


= n[x 


r, t] /. {coeff[[8]]}; 


ort4 [xO 


x_] 


= - 1 / (D [po4 [x] , x] ) * (x - xO) + po7 [xO] ; 



*) 



ort7[x0_, x_] = -1/ (D[po7[x], x] ) * (x - xO) +po7[x0]; 
ortS [xO_, x_] = - 1 / (D [po8 [x] , x] ) * (x - xO) + poS [xO] ; 



ln[18]:= Table [ 

intersectortSpoV [m] = 

FindInstance[{Evaluate[ort8[0.25 + 0.1 *m, x] ] ==po7[x], x>0.18, x<0.8}, x. Reals] [ [1] ][ [ 
1]] (* determine intersections *) 

'{m, 
0, 
4}] 



Out[18]= {x ^ 0.2382, x -> 0.3202 95, x ^ 0.403032, x ^ 0.487458, x ^ 0.574101} 
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ln[19]:= Table [ 

intersectort7po8old[n] = 
Findlnstance [ {Evaluate [ort7 [0.25 + 0. 1 * n, x] ] == po8 [x] , x > 0.18, x < 0.8}, x. Reals] [ [1] ] [ [ 
1]] 

'{n, 
0, 
4}] 

Out[19]= {x^O. 279988, 0.40 9555, x^O. 530964, x^O. 645, x^O. 754393} 

ln[20]:= Show[Plot [ { (*neu8 [x] , neu7 [x] , *) po7 [x] , po6[x], po5[x], po4[x], 

po3[x] , po2[x] , pol[x] , Abs[x] , 0}, {x, 0, 0.8}, PlotRange -» {-0.05, 1.}, 
AspectRatio -» Automatic, Background -> White, PlotStyle -» {{Gray, Thick}}, Axes -» False] , 
Table[Plot [Evaluate [ort7 [0.2 + 0.09 * m, x] ] , {x, 0.17, 0.8), 

AspectRatio -» Automatic, Background -» White, PlotRange-* {0., 1.}, 
PlotStyle^ {Black, Dashed, Thick}, Axes -» False] , {m, 0, 8}] 
(*Table [Plot [Evaluate [y8 [x/ . Schnitty7mitneu8 [n],x]],{x,0.2,0.8} , AspectRatio->Automatic, 
Background-»White,PlotRange-»{0 .4, 1 .2} ,PlotStyle-»{Blue} , Axes-»False] , {n, 0,3}],*) 



Out[20]= 




■ Normal Vectors 

ln[21]:= slope? [x_] = -1/ (D[po7[x], x] ) 

Out[21l= |{- 

^ 0.282592 - 1.72044 x- 0.547624 x^ 




In[22]:= 



vecend[n_] = {0.25 + 0.1*n, po7 [0 . 25 + . 1 * n] [[1]] [[1]]} + 

{0.25 + 0.1*n + l, po7[0.25 + 0.1*n] [[1]] [[1]] + slope? [0 . 25 + . 1 * n] [[1]] [[1]]} 



Out[22]= 



{ 




0.565184 (0.25 + 0. In) - 1.72044 (0 . 25 + . 1 n) ^ - . 365083 (0.25 + O.ln)^ 



ln[23]:= 



slopes [x_] = -1/ (D[po8[x], x] ) 



{{ 



1 



ill 



Out[23]= 



0.287713-1. 4995 6x-0. 424976 X- 



ln[24]:= Table [ 

intersectort7po8 [n] = x / . Findlnstance [ 

{Evaluate [ort7[ 0.25 + 0.1 *n, x] ] == po8 [x] , x > 0.18, x < 0.8}, x. Reals] [ [1] ] [ [1] ] 



4}] 

Out[24]= {0.279988, 0.409555, 0.530964, 0.645, 0.754393} 

ln[25]:= Table [vecendS [m] = 

{intersectort7po8 [m] , po8 [intersectortVpoS [m] ] [ [ 1] ] [ [ 1] ] } + { intersectortVpoS [m] + 1, 
po8[intersectort7po8[m] ] [ [1] ] [ [1] ] + slope8 [intersectort7po8 [m] ] [ [ 1] ] [ [ 1] ] } , {m, 0, 

Out[25]= {{1.55998, 7. 9665}, {1.81911, 4 .36444}, {2.06193, 3.31738}, {2.29, 2 .72453}} 



{n, 
0, 
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ln[26]:= Show [Plot [ {po8 [x] , po7 [x] , po6[x], po5 [x] , ftbs [x] } , {x, 0.1, 0.9}, PlotRange -> {0.5, 1.3}, 
AspectRatlo -> Automatic, Background -» White, PlotStyle -» {{Gray, Thick}}, Axes -* False] , 
Table [Plot [Evaluate [ort7[0. 25 + 0.1 *m, x] ] , {x, 0.2, 0.9}, 

AspectRatio -» Automatic, Background -» White, PlotRange -» {0.5, 1.2}, 
PlotStyle^ {Gray, Dashed, Thick}, Axes ^ False] , {m, 0, 3}], 
Table [Plot [Evaluate [ortS [x /. intersectort7po8old[n] , x] ] , {x, 0.2, 0.8}, 
AspectRatio -» Automatic, Background -» White, PlotRange ^ {0. 5, 1.2}, 
PlotStyle-* {Gray, Dashed, Thick}, Axes ^ False] , {n, 0, 3}], 
Table [Graphics [{Thick, Black, Arrow [{{0.25 + 0.1 *n, po7[0.25 + . 1 * n] [ [1] ] [ [1] ] } , 
{0.25 + 0.1*n + 0.7 / Norm[vecend[n] ] ^2 * 1, po7[0.25 + 0.1* n] [[1]] [[1]] + 
0.7 / Norm[vecend[n] ] ^2 * slope7 [0 . 25 + 0.1 * n] [ [ 1] ] [ [ 1] ] } } ] } ] , {n, 0, 3}] , 
Table [Graphics [ {Thick, Black, Arrow [{ {intersectort7po8 [n] , 
po8 [intersectort7po8 [n] ] [ [1] ] [ [1] ] } , 
{intersectort7po8 [n] + . 7 / Norm [vecend8 [n] ] ^ 2 * 1, po8 [intersectort7po8 [n] ] [ [1] ] [ [1] ] + 
0.7 / Norm[vecend8[n] ] *2 * slope8[intersectort7po8[n] ][ [1] ][ [1] ]}}]}] , {n, 0, 3}] 



Out[26]= 




Polar Non-Orthogonal Grid Ansatz 

Generation of Quasi-Polar Grid in 2D 



Ansatz for Quasi-Polar Coordinates 



■ Coordinate Transformation with Product Ansatz 



We define a coordinate transformation with product ansatz, which is reasonable for any kind of coordinate system. 

In[1]:= X = X[f , ?7] ; 

y = n]} 
ln[3]:= X[€'_, J7_] = a[f] * a[ri] ; 

ln[5]:= a [77_] = (1 + alphal * 77 + alpha2 * tj * 2 + alphaS * 77 ^ 3) * Cos [77] ; 

= (1 + betal * 77 + beta2 * r7''2 + betas * J7'^3) * Sin[j7] ; 
a[f_] = (al* f + a2*f*2) ; 
b[f_] = (bl* S*b2*S^2) ; 

We calculate the Jacobian of this transformation to get the base vectors in the new coordinate system. 

In[9]:= lij = {{D[x, f] , D[x, rj] } , {D[y, f ] , D[y, 77]}}; 
ln[iO]:= li jtransposed := Transpose [lij] 



■ Covariant Components of lUletric Tensor 

For trivial signatures, the metric tensor is simply given by the Jacobian times its transpositon. The off diagonal elements contain 
the information about the non-orthogonality (angles), the diagonal elements are measures of lenght in the coordinate directions. 

In[11]:= gij = li jtransposed. lij; 

ln[12l:= gll[f_, J7_] =gij[[i. 1]]; 
gl2[f_, r,_] = gij[[l, 2]]; 

g2i[f_, 77_] = gij [[2, 1]]; 

g22[f_, J7_] = gij [[2, 2]]; 



■ Covariant Base Vectors 



The base vectors yield 

ln[16]:= ef[f_, J7_] = lijtransposed[ [1] ] 

Out[16]= { (1 + alphal r] + alpha2 rj^ + alpha3 rj^) (al + 2 a2 .f) Cos [rj] , 
(1 +betal rj +beta2 17^ +beta3 17^) (bl + 2 b2 ^) Sin [17]} 
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ln[17]:= &n[€-, 1—] = lijtransposed[ [2] ] 

Out[17]= I (alphal + 2 alpha2 r; + 3 alphas r;^ ) (al + a2 cj^ ) Cos [ 17] - 

(1 + alphal r7 + alpha2 r]^ + alphas r]^) [al S + sl2 S^) Sin[f7], 
(1 + betal T] +beta2 rj^ + betaS rj^) (bl S + b2 f^) Cos [rj] + 
(betal + 2 beta2 17 + 3 beta3 rj^) (bl S + h2 S^) Sin [f]] } 



PDEs + Boundary Conditions 



■ Controling the Grid 

In order to control the shape of the grid by a manageable number of parameters and for reasons of symmetry we impose geometri- 
cally motivated boundary conditions. Since some of them are either redundant or exclusive, the following list is somewhat 
heuristic. 

In[18]:= ESSol = Solve [ { 

(*gl2 [f , rj] ==0, *) {* orthogonality *) 

e{[e, 0][[1]] ==1, (* scaling *) 

ef[f, 0][[2]] ==0, {* in "x"-axis *) 

Gf [f. Pi / 2] [ [1] ] ==0, {* in "y"-axis *) 

eslS, Pi/ 2] [[2]] ==1, (* scaling *) 

^rilSi 0] [[!]] == 0, (* normal to "x"-axis *) 

ejf, 0] [[2]] ---- S, {* scaling *) 

[f , Pi / 2] [ [1] ] == -f, {* orientation *) 
e^ilS, Pi/ 2] [[2]] ==0, (* normal to "y"-axis *) 
X[0, 0] == 0, (* set origin *) 
Y[0, 0] == 0, {* set origin *) 
X[f, 0] == X[f, 2* Pi], {* periodicity *) 
Y[f, 0] == Y[f, 2* Pi], {* periodicity *) 
X[f, 0] == f*X[l, 0] (* scaling *) 

}' 

{alphal, alpha2, betal, beta2, , alphaS, betaS, al, a2}] 

Solve ::svars : Equations may not give solutions for all "solve" variables. » 

beta3 71^ - 1 6 b2 f + b2 beta3 71^ S 
Out[18]= {{betal ^ , 

4 7T (1 + b2 f ) 

-betas jr^ + 4 b2 f - b2 betaS /r^ S 

beta2 ^ , alpha2 ^ 0, alphal ^ 0, alphaS -> 0, al ^ 1, a2 ^ O}} 

(l+b2 S) 

ln[19]:= Xnew[5'_, ri_] = Evaluate [X[f, 77] /. ESSol[[l]]] //Simplify 
Ynew[f_, ri_] = Evaluate [Y [f, 77] /. ESSol[[l]]] // Sin^jlify 

Out[19]= fCos[r7] 
1 

Out[20]= f (bl + b2 f ) 

4 TT^ (1 +b2 S) 

(- 16 b2 77 r] f + 16 b2 rj^ f + betaS 7:" rj ( 1 + b2 f ) - 4 betaS 7T^ rj^ (1 + b2 + 4 71^ (l + betaS rj^) (1 + b2 f ) ) 
Sin[)7] 



■ Plots 



This example for a quasi-polar mesh is being controled by three parameters. 



quasi_polar_grid.nb 1 3 



ln[21]:= Manipulate [ 

ParametricPlot [-[f Cos [rj] , 



Out[21]= 



4 TT^ (1 +b2 S) 



f (bl + b2 f ) (- 16 b2 TT 77 f + 16 b2 77^ f + betaS tt^ rj (1 + b2 f ) 



4beta3 7r^r7^ (l + b2f) + 4 tt^ (l+betaSrj^) (l+b2f)) Sin [77]}, 
{S, 0, 1}, {77, 0, Pi / 2}, AspectRatio -» 1, PlotRange -» {0, 1}, Axes ^ None, 
BoundaryStyle -» {RGBColor [ .2, .2, .2], Thick}, 

MeshStyle -» { {RGBColor [. 2 , .2, .2], Thick}}, MeshFunctions -» Automatic, 
MeshShading -> {{White, White}, {White, White}}, Background -» White (*, LabelStyle-**) ] , 
{{bl, 1}, -1, 1}, {{b2, 0}, -1, 1}, { {betas, 0}, -1, 1}] 




Geometrical Characteristics of Toy Models 

Determination of maximal relative volumes respectively off-diagonal metiic elements. 

In[221:= li jnew = { {D [Xnew [ f , r?] , f] , D [Xnew[f , 77] , 77]}, {D [Ynew [f , 7?] , S] , D [Ynew , 77] , 77]}}; 
ln[23]:= gijnew = Transpose [lijnew] .lijnew; 



ln[24]:= grida = {bl -» 1, b2 -» 0, betaS -* 0}; 

gridb = {bl ^0.8, b2 -» 0, beta3 -» 0} ; 
gride = {bl -» 0.8, b2 -» 0, beta3 -» 0.5}; 
gridd = {bl -» . 8, b2 -» . 2, betaS -» 0.5}; 



4 I quasi_polar_grid.nb 



ln[28]:= deta[€_, ri_] = Evaluate [Det [gijnew] /. grida] ; 

detb[€'_, 77_] = Evaluate [Det [gijnew] /. gridb] ; 

detc[f_, r]_] = Evaluate [Det [gijnew] /. gride]; 

detd[f_, ri_] = Evaluate [Det [gijnew] /. gridd] ; 

ln[32]:= offa[€_, 77_] = Evaluate[gijnew[ [1, 2]] /. grida]; 

offb[€_, 7j_] = Evaluate[gljnew[ [1, 2]] /. gridb]; 

offc[S_, 77_] = Evaluate[gijnew[ [1, 2]] / . gride]; 

offd[€_, 77_] = Evaluate[gijnew[ [1, 2]] /. gridd]; 

■ Maxima 

ln[36]:= FindMaxim™ [ {deta [f , rj] / f 2, fa S& 77 a && f s 1 &S ^7 s Pi / 4}, {f, r)}] 

Out[36]= {1., {f^ 0.899989, 17^ 0.706847 }} 

ln[37]:= FindMaximum [ {Abs [ of f a [f , rj] / S] , fa 0&S77£OSSfSlS&?7S Pi/4}, {f, 77}] 

Out[371= {0., {S ^ 0., r] ^ 0.]] 

In[38]:= FindMaximum [ {detb [ f , rj] / f ^2, fa SS ^ £ S& f S .8 SS s Pi / 4}, {f, 77}] 

Out[38]= {0.64, {f^ 0.719994, r] ^ . 706852} } 

ln[39]:= FindMaximvim [ {Abs [of fb [ f , 77] / f] , f £ SS 77 a SS f S 1 SS 77 S Pi / 4}, {f, 77}] 

Out[39]= {0.179999, { 4" . 85142 9, j] ^ . 78387 9} } 

ln[40]:= FindMaxiinum[{detc[f, 77] / f *2, fa SS 77 £ SS f S 1 SS 77 S Pi / 4}, {f , 77}] 

Out[40]= {1.19161, {f^ 0.860374, 77^ 0.327134}} 

ln[41]:= FindMaximum[{Abs[offc[f, 77] / f ] , fa SS 77 a SS f S 1 SS 77 S Pi / 4}, {f, 77}] 

Out[41l= {0.128795, {f^ 0.899702, 77^ 0.785398}} 

ln[42]:= FindMaxiinum[{detd[f, 77] / f *2, fa SS 77 a SS f S 1 SS 77 S Pi / 4}, {f, 77}] 

Out[42]= {1. 62084, {f ^ 1., 77 ^ 0.293574}} 

ln[43]:= FindMaximum[{Abs[offd[f, 77] / f ] , fa SS 77 a SS f i 1 SS 77 S Pi / 4}, {f, 77}] 

Out[43]= {0.248654, {f^ . 999999, 77 ^ 0.424638}} 

■ Minima 

ln[44]:= FindMinimum[{deta[f, 77] / f *2, fa SS 77 a SS f S 1 SS 77 S Pi / 4}, {f, 77}] 

Out[44]= {1., {f^ . 899989, 77 ^ 0.706847}} 

ln[45]:= FindMinimum[{Abs[offa[f, 77] / f ] , fa SS 77 a SS f S 1 SS 77 S Pi / 4}, {f, 77}] 

Out[45]= {0., {f^O., 77^0.}} 

In[46]:= FindMinimum[{detb[f , 77]/f*2, fa 0&&77aO&&f£ .8&&77:£ Pi/4}, {f, 77}] 

Out[461= {0.64, {f -> 0.719994, 77 ^ 0.706852}} 



In[47]:= FindMinimum[{Abs[offb[f, rj] / S] , S ^ SS i SS f i 1 SS S Pi / 4}, {S, rj}] 

Out[47]= {8.29014 X 10"\ {f^O. 772027, r] ^ 2 . 30282 x lO"' } } 

ln[48]:= FindMinimum[{detc[f, 77] / f*2, fa SS 17 £ SS f S 1 SS 77 S Pi / 4}, {f, I?}] 

Out[48]= {0.757631, {f^ 0.897993, ?7 ^ . 785396} } 

ln[49]:= FindMinimum[{Abs[offc[f, ??] / f ] , S ^ SS rj i SS f S 1 SS 77 S Pi / 4}, rj}] 

Out[49]= {2 .46448 X 10"\ {f ^ . 734207, r] ^ 6 . 84577 x 10"*} } 

ln[50]:= FindMinimum[{detd[f, 77] / f 2, fa SS 77 £ SS f S 1 SS 77 S Pi / 4}, {f , 77}] 

Out[50]= {0.75763, {f ^ 5 . 89009 x 10"% r] ^ . 785398} } 

ln[53]:= FindMinimum [ {Abs [of f d [f , 77] / f ] , S h && 77 a && f :£ 1 && 77 £ Pi / 4} , {f , 77}] 

FindMinimum : : eit : 
The algorithm does not converge to the tolerance of 4 . 806217383937354 ~ 
in 500 iterations. The best estimated solution, with 
feasibility residual, KKT residual, or complementary residual 
of {0.000658608, 0.100761, 0.000246563}, is returned. » 

Out[53]= {3.42416 X 10"S {f^ 0.709494, r] ^ . 716048} } 
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Abstract 



We study conservative numerical methods for implicit radiation hydrodynamics (RHD) 
in general non-steady coordinates in 2D and 3D. 

We discuss fundamental mathematics and physics of RHD with special focus on dif- 
ferential geometric consistency and stTidy numerical methods for nonlinear conservation 
laws. We apply Vinokurs theorem to obtain the strong conservation form for conser- 
vation laws in general curvilinear coordinates. Differential geometric derivations in this 
context lead to a reformulation of artificial viscosity in such general coordinates. 

A non-static gravitation equation in order to avoid solving the Poisson equation for 
the gravitational potential is suggested. 

It is shown that there is no coordinate system in 2D and 3D that contains the po- 
lar respectively spherical coordinates and allows multi-dimensionally adaptiveness. We 
suggest an ansatz to generate non-orthogonal non-steady coordinates that contain these 
orthogonal reference frames. 

Zusammenfassung 

Wir untersuchen konservative, implizite numerische Methoden im Rahmen der Strahlung- 
shydrodynamik (SHD) in allgemeinen Koordinaten. 

Grundlegende mathematische und physikalische Konzepte werden mit speziellem Fokus 
auf geometrische Konsistenz diskutiert und numerische Methoden fiir nichtlineare Er- 
haltungssatzc untcrsucht. Mit Hilfe des Vinokurschen Theorems lasst sich die konser- 
vative Form der Erhaltungssatze auf allgemeinen Koordinaten formulieren. Differential- 
geometrische Uberlegungen in diesem Zusammenhang fiiliren zu einer Umformulierung 
der kiinstlichen Viskositat fiir solch allgemeine Koordinatensysteme. 

Wir schlagen eine nicht-statische Gleichung fiir die Eigcngravitation vor, um die L6- 
sung der Poisson-Gleichung fiir das Gravitationspotential zu umgehen. 

Es wird gezeigt, dass kein Koordinatensystem in 2D und 3D existiert, das die Polar- 
bzw. Kugelkoordinaten beinhaltet und mehrdimensionale Adaptivitat des Gitters er- 
laubt. Wir schlagen einen Ansatz fiir nicht-orthogonale zeitabhangige Koordinaten vor, 
der diese Referenzgitter enthalt. 



Curriculum Vitae 



Name: 

Datum und Ort der Geburt: 
Eltern: 



Harald HOLLER Bakk.rer.nat 

04. Marz 1983 in Villach, Osterreich 

Anna HOLLER und Engelbert Viktor HOLLER 



Adresse: 

E-Mail: 

Homepage: 



Alszeile 119, 1170 Wien 

harald . hoellerSunivie . ac . at 

http : / /homepage . univie . ac . at/harald . hoeller/ 



Stationen 



2009, 2010 
08/2007-01/2010 



2007, 2008 



2007, 2008, 2009 



Sommer 2005, 2006 



Lektor an der Fakultat fiir Physik - Lehrveranstaltungen aus The- 
oretischer Physik und Mathematischer Methoden. 
Projektmitarbeiter bei eLearnPhysik - E-Learning an der Fakultat 
fiir Physik. Schwerpunkte waren die redaktionelle und technische 
Betreuung des Wiki der Fakultat fiir Physik, Content und Usabil- 
ity Fragen, sowie E-Lcarning im Bereich Theoretische Physik und 
der Einsatz von Computeralgebra in der Lehre. 
eTutor an der Fakultat fiir Physik - Lehrveranstaltungen aus Ex- 
perimentalphysik sowie Theoretischer Physik und Mathematischer 
Methoden. 

Seminare und Workshops im Rahmcn der KindcrUni Wien, 
des Wissenschaftsclub fiir Jugendliche sowie der StauneLaune 
Forschungswochen. 

Ferialpraktikant bei Infineon Technologies Austria AG am Stan- 
dort Villach im Bereich Forschung und Entwicklung. 



Schule und Studium 



seit 2005 Magisterstudium Astronomie. 

seit 2003 Diplomstudium Physik (Zweitstudium) an der Universitat Wien. 
2002, 2005 Zuerkennung von Leistungsstipendium. 

2005 Abschluss des Bakkalaureatsstudiums Astronomie mit ausgezeich- 
netem Erfolg. 

2005 Unterstellung in neuen Studienplan; Bakkalaureatsstudium As- 
tronomie. 

2001-2005 Diplomstudium Astronomie an der Universitat Wien. 

2001 Matura am BORG Klagenfurt mit ausgezeichnetem Erfolg. 



Publikation 



• Wiki Based Teaching and Learning Scenarios at the University of Vi- 
enna, Harald Holler und Peter Reisinger, Vortrag auf der World Conference on 
Educational Multimedia, Hypermedia & Telecommunications (ED-MEDIA 30.06.- 
04.07.08, Technische Universitat Wien), in den Proceedings erschienen. 



Eingeladene Vortrage 

• Hearing zum Finale des Medida Prix 2009, Franz Embacher, Harald Holler, 
Christian Primetshofer, Peter Reisinger, Brigitte Wolny - Vorstellung des Projekts 
eLeamPhysik, 15.09.2009, Berlin. 

• Phaidra, Wiki und die Content-Strategie der Fakultat fiir Physik, Harald 
Holler, Peter Reisinger im Rahmen des Phaidra Day an der Fakultat fiir Physik 
am 01.10.2009 . 

• eLearnPhysik, Franz Embacher, Harald Holler, Christian Primetshofer, Peter 
Reisinger - Vorstellung des Projekts im Rahmen der Friday Lectures am Projek- 
tzentrum Lehrentwicklung am 07.12.2007. 



